44

Is there an easy way of transforming Shapely objects (namely, Polygons and MultiPolygons) from one projection to another without having to dig around and extract coordinates by hand?

In fact, I don't even care if they are Shapely objects at this point, I just want to pass features and a projection, and get a reprojected set of features back.

Does this sort of functionality exist, or must it be hand coded?

1
  • 2
    I believe that's out of the scope of Shapely, you may want to look into Fiona. fiona.transform looks to have what you need. Dec 22, 2014 at 21:51

3 Answers 3

71

While shapely doesn't natively understand coordinate systems, shapely.ops.transform() can do that along with pyproj. If pyproj.Proj can understand your both of your coordinate systems, then it can be made into a function that shapely can transform with.

From the shapely docs:

pyproj >= 2.1.0

import pyproj

from shapely.geometry import Point
from shapely.ops import transform

wgs84_pt = Point(-72.2495, 43.886)

wgs84 = pyproj.CRS('EPSG:4326')
utm = pyproj.CRS('EPSG:32618')

project = pyproj.Transformer.from_crs(wgs84, utm, always_xy=True).transform
utm_point = transform(project, wgs84_pt)

pyproj < 2.1

from functools import partial
import pyproj
from shapely.ops import transform

project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:4326'), # source coordinate system
    pyproj.Proj(init='epsg:26913')) # destination coordinate system

g2 = transform(project, g1)  # apply projection
4
  • 6
    If you do not want to use the itertools module you can do project = lambda x, y: pyproj.transform(pyproj.Proj(init='epsg:4326'), pyproj.Proj(init='epsg:26913'), x, y) and then g2 = transform(project, g1).
    – Elmex80s
    Jan 24, 2018 at 13:15
  • 2
    This suggested answer is for pyproj1 whereas there is now a preferred transformation for pyproj2 using Transformer. See here: pyproj4.github.io/pyproj/stable/gotchas.html Jul 9, 2019 at 14:59
  • 1
    I found Nick O'Dell's answer below (gis.stackexchange.com/a/328642) to be much faster. Dec 27, 2019 at 4:21
  • Checking my knowledge: wgs84_pt describes the coordinates 72.2495W (longitude) and 43.886N due to the always_xy=True. For always_xy=False lat/lon would need to be swapped when constructing the initial point. Is this correct?
    – Dave
    Mar 8, 2021 at 20:47
26

If you're using pyproj2, it's much easier to use a Transformer. Here's an example:

import pyproj
from shapely.ops import transform

project = pyproj.Transformer.from_proj(
    pyproj.Proj(init='epsg:4326'), # source coordinate system
    pyproj.Proj(init='epsg:26913')) # destination coordinate system

# g1 is a shapley Polygon

g2 = transform(project.transform, g1)  # apply projection

This is also much faster, becase pyproj does not need to recreate the projection for every point.

2
  • 3
    Way faster than the accepted answer! Dec 27, 2019 at 4:20
  • Really much faster! fantastic.Great!
    – Jane
    Jan 9, 2020 at 13:12
14

While not a Shapely solution, using GeoPandas allows for relatively straightforward projection. For example, if we want to convert a shapefile to ESPG 4326:

import geopandas as gpd

HabModelEnviro = gpd.GeoDataFrame.from_file('data/HabModelEnviro.shp').replace({-999: None})

HabModelEnviroWGS84 = HabModelEnviro.to_crs({'proj':'longlat', 'ellps':'WGS84', 'datum':'WGS84'})
1
  • 8
    Geopandas uses Shapely (look at geodataframe.py for example)
    – gene
    Dec 23, 2014 at 16:34

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.