1

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.

  • 1
    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. – Vince Dec 30 '19 at 12:39
  • 1
    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 – Zoltan Dec 30 '19 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. – jelleve Dec 30 '19 at 13:00
1

@jelleve, your answer is correct and in line with https://pyproj4.github.io/pyproj/stable/advanced_examples.html#optimize-transformations.

Additionally, you may be interested in geopandas. It will enable you to read/write these geometries and reproject them using the same method you showed above using one line of code:

gdf = geopandas.read_file("input.geojson")
gdf_wgs84 = gdf.to_crs("epsg:4326")
gdf_wgs84.to_file("output.geojson", driver="GeoJSON")

See also: https://geopandas.readthedocs.io/en/latest/projections.html

0

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.

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