Project large collection of shapely object to different CRS [closed]

What is the most efficient way of projecting several thousands of polygons and multipolygons to a different CRS in Python? I have found that looping over the polygons and reprojecting them one by one is rather slow. Ideally, for a thousand polygons, this operation shouldn't take longer than 1 second.

Currently, I am using an approach where I extract all coordinates into a single numpy array, convert them, and split the results back up into polygons. While this is quite fast, right now it doesn't handle polygons with interiors or multipolygons. I could extend this approach, but before doing so, I was wondering whether I am overlooking something that is a bit more manageable.

• You can have it good, cheap, or fast (pick any two). Unfortunately, this question is more of an open-ended discussion topic, more suited to Geographic Information Systems Chat than the Focused question/Best answer model used in our main site. It's difficult to tell if you want a software recommendation, algorithmic help, or something else. It's also unclear how many vertices your features contain, and what the exact CRS are involved, to gauge whether your goal is even feasible. Please Edit the Question to add details and a specific question. Commented Dec 30, 2019 at 12:39
• What about ogr2ogr? You can call it from Python or use ogr and osr from osgeo Python package. see: gis.stackexchange.com/questions/154004/… or pcjericks.github.io/py-gdalogr-cookbook/projection.html Commented Dec 30, 2019 at 12:39
• Thank you for your inputs. I have managed to solve the problem, which was nothing more than an error on my part. Commented Dec 30, 2019 at 13:00

After some experimenting, it seems that I was just plain wrong. My initial code looked as follows:

def convertCoordinates(polygon):
transformer = pyproj.Transformer.from_proj(pyproj.Proj(init='epsg:27563'), pyproj.Proj(init='epsg:4326'))
projected = shapely.ops.transform(transformer.transform, polygon)
return projected

for polygon in polygons:
projected_coordinates = convertCoordinates(polygon)


As this was quite slow, I concluded that the loop over the polygons must be the bottleneck. This was confirmed after implementing my second approach (described in the original question), which was much faster.

However, it turns out that the loop itself was not the bottleneck. The coordinate transformation seems to be a reasonably light-weight operation. The creation of the Transformer object, however, is not. As such, it is a bad idea to define and initialize it multiple times. The solution below seems to result in the right combination of speed and elegance:

def convertCoordinates(polygon, transformer):
projected = shapely.ops.transform(transformer.transform, polygon)
return projected

transformer = pyproj.Transformer.from_proj(pyproj.Proj(init='epsg:27563'), pyproj.Proj(init='epsg:4326'))
for polygon in polygons:
projected_coordinates = convertCoordinates(polygon, transformer)


This version handles +- 700 polygons of various complexity in +-0.2 seconds.

gdf = geopandas.read_file("input.geojson")