I would like to project all the features in a geojson from WGS84 (EPSG:4326) to Web Mercator (EPSG:3857) using python.

This geojson contains only polygons.

I found out the following shapely code combined with pyproj will do it. However after the projection has been applied then I need to manually build the output geojson with the new Web Mercator projection.

import pyproj
import json

from functools import partial

# on windows: conda install -c conda-forge shapely=1.5.17
from shapely.geometry import Polygon, Point, shape, mapping
from shapely.ops import transform

# the shapely projection function
datum_wgs84 = pyproj.Proj(init='EPSG:4326') 
projection_web_mercator = pyproj.Proj(init='EPSG:3857')
projection_wm_func = partial(pyproj.transform, datum_wgs84, projection_web_mercator)

# load GeoJSON file containing sectors
with open('my_input_file.geojson') as f:
  js = json.load(f)

# project GPS lat/lon coordinates to web mercator for each polygon in the geojson
for feature in js['features']:
  polygon = shape(feature['geometry'])
  print('{} bounding box (in WGS84): {}'.format(feature['properties']['my_polygon_name'], polygon.bounds))
  print('-> A sample pair in the original WGS84 projection: {}'.format(mapping(polygon)['coordinates'][0][0]))
  projected_polygon = transform(projection_wm_func, polygon)
  print('-> A sample pair in the projected "Web Mercator": {}'.format(mapping(projected_polygon)['coordinates'][0][0]))

# TODO put all the projected polygons together in a new output geojson

Is there a more compact way of doing this manual part to avoid most of the boiler plate code? Ideally some sort of function with inputs:

  • input geojson filename
  • output geojson filename
  • existing projection
  • final projection

Or at least a combination of functions to avoid the for loop.

  • What's wrong with 'for'? – nickves Jul 19 '17 at 14:55

This can be done without shapely. See below, adapted from this example.

from pyproj import Proj
import json
from pprint import pprint

in_file = '/path/to/fc.geojson'

p_web = Proj(init='EPSG:3857')

with open(in_file) as src:
    fc_in = json.load(src)

# Define dictionary representation of output feature collection
fc_out = {'features': [],
          'type': 'FeatureCollection'}

# Iterate through each feature of the feature collection
for feature in fc_in['features']:
    feature_out = feature.copy()
    new_coords = []
    # Project/transform coordinate pairs of each ring
    # (iteration required in case geometry type is MultiPolygon, or there are holes)
    for ring in feature['geometry']['coordinates']:
        x2, y2 = p_web(*zip(*ring))
        new_coords.append(zip(x2, y2))
    # Append transformed coordinates to output feature
    feature_out['geometry']['coordinates'] = new_coords
    # Append feature to output featureCollection


You could try using ogr2ogr utility:

ogr2ogr -f "GeoJSON" output.geojson input.geojson -s_srs EPSG:4326 -t_srs EPSG:3857 

This will not be the final answer, but i will try to help you:

1) If you want a compact way just put the code in a python function. foor loops always are needed.

2) In order to put all the projected polygons together in a new output geojson, only must to build a dictionary and write wherever you want. in the next code i use ogr library:

    print ("Begin")

features =[]

# loop through the features in the layer
feature = layer.GetNextFeature()
while feature:
  counter +=1
  if counter > 10 : break


  # destroy the feature and get a new one
  feature = layer.GetNextFeature()

geo_json_dict = {"type": "FeatureCollection",
                      "features": features

# open the output text file for writing
file = open('features.json', 'w')

# close the data source and text file
print ("End")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.