Software developer

How to convert geojson from polar to mercator projection

Structure of geojson file is quite simple. To work with geojson we can use standart module in python called json. Often geojson files contains utf-8 symbols, so codecs module will be helpfull. And main module that we need for projections is pyproj (follow installation guide or simply run command below).

pip install pyproj

Pyproj provides APIs to convert pair of coordinates. Projection of single point will look like:

polar = pyproj.Proj(init='epsg:4326')
mercator = pyproj.Proj(init='epsg:3857')

[x, y] = [26.2513165, 50.6196175]
new_x, new_y = pyproj.transform(polar, mercator, x, y)

All we need to convert any geomentry in geojson is simple recursive function:

def project_coords(coords, from_proj, to_proj):
    if len(coords) < 1:
        return []

    if isinstance(coords[0], numbers.Number):
        from_x, from_y = coords
        to_x, to_y = pyproj.transform(from_proj, to_proj, from_x, from_y)
        return [to_x, to_y]

    new_coords = []
    for coord in coords:
        new_coords.append(project_coords(coord, from_proj, to_proj))
    return new_coords

Actually you can use this example to convert geojson files with different projections. More details about this you can find in pyproj documentation.

Full code snippet of convertation from Polar projection (EPSG:4326) to Pseudo-Mercator (WGS84 or EPSG:3857) and vise versa:

import json
import codecs
import pyproj
import numbers

# Define projections
polar = pyproj.Proj(init='epsg:4326')
mercator = pyproj.Proj(init='epsg:3857')


def project_coords(coords, from_proj, to_proj):
    if len(coords) < 1:
        return []

    if isinstance(coords[0], numbers.Number):
        from_x, from_y = coords
        to_x, to_y = pyproj.transform(from_proj, to_proj, from_x, from_y)
        return [to_x, to_y]

    new_coords = []
    for coord in coords:
        new_coords.append(project_coords(coord, from_proj, to_proj))
    return new_coords


def project_feature(feature, from_proj, to_proj):
    if not 'geometry' in feature or not 'coordinates' in feature['geometry']:
        print('Failed project feature', feature)
        return None

    new_coordinates = project_coords(feature['geometry']['coordinates'], from_proj, to_proj)
    feature['geometry']['coordinates'] = new_coordinates
    return feature


def read_data(geom_file):
    with open(geom_file, encoding='utf-8') as data:
        data = json.load(data)
    return data


def save_data(data, geom_file):
    json_data = json.dumps(data, indent=2)
    f = codecs.open(geom_file, "w")
    f.write(json_data)
    f.close()


def project_geojson_file(in_file, out_file, from_proj, to_proj):
    data = read_data(in_file)

    projected_features = []
    for feature in data['features']:
        projected_features.append(project_feature(feature, from_proj, to_proj))
    data['features'] = projected_features

    save_data(data, out_file)


def project_file_to_polar(in_file, out_file):
    project_geojson_file(in_file, out_file, from_proj=mercator, to_proj=polar)


def project_file_to_mercator(in_file, out_file):
    project_geojson_file(in_file, out_file, from_proj=polar, to_proj=mercator)


# Project geojson file from mercator to polar
project_file_to_polar('example-mercator.geojson', 'example-polar.geojson')

# And from polar to mercator
project_file_to_mercator('example-polar.geojson', 'example-mercator.geojson')

# Or any other projection you want
project_geojson_file('example-epsg-6069.geojson', 'example-polar.geojson',
                     from_proj=pyproj.Proj(init='epsg:6069'), to_proj=polar)

* All code in this article provided under MIT license

Thank you for attention!