added ability to project geometries, geodataframes, and graphs to any…

… CRS

osmnx/projection.py

@@ -15,7 +15,7 @@
 from .utils import log, make_str
 
 
-def project_geometry(geometry, crs, to_latlong=False):
+def project_geometry(geometry, crs, to_crs=None, to_latlong=False):
     """
     Project a shapely Polygon or MultiPolygon from lat-long to UTM, or vice-versa
     
@@ -25,8 +25,10 @@ def project_geometry(geometry, crs, to_latlong=False):
         the geometry to project
     crs : dict
         the starting coordinate reference system of the passed-in geometry
-    to_latlong : if True
-        project from crs to lat-long, if False, project from crs to local UTM zone
+    to_crs : dict
+        if not None, just project to this CRS instead of to UTM
+    to_latlong : bool
+        if True, project from crs to lat-long, if False, project from crs to local UTM zone
     
     Returns
     -------
@@ -38,12 +40,12 @@ def project_geometry(geometry, crs, to_latlong=False):
     gdf.gdf_name = 'geometry to project'
     gdf['geometry'] = None
     gdf.loc[0, 'geometry'] = geometry
-    gdf_proj = project_gdf(gdf, to_latlong=to_latlong)
+    gdf_proj = project_gdf(gdf, to_crs=to_crs, to_latlong=to_latlong)
     geometry_proj = gdf_proj['geometry'].iloc[0]
     return geometry_proj, gdf_proj.crs
 
 
-def project_gdf(gdf, to_latlong=False):
+def project_gdf(gdf, to_crs=None, to_latlong=False):
     """
     Project a GeoDataFrame to the UTM zone appropriate for its geometries' centroid. The simple calculation
     in this function works well for most latitudes, but won't work for some far northern locations like
@@ -52,7 +54,9 @@ def project_gdf(gdf, to_latlong=False):
     Parameters
     ----------
     gdf : GeoDataFrame
-        the gdf to be projected to UTM
+        the gdf to be projected
+    to_crs : dict
+        if not None, just project to this CRS instead of to UTM
     to_latlong : bool
         if True, projects to latlong instead of to UTM
     
@@ -63,48 +67,56 @@ def project_gdf(gdf, to_latlong=False):
     assert len(gdf) > 0, 'You cannot project an empty GeoDataFrame.'
     start_time = time.time()
 
-    if to_latlong:
-        # if to_latlong is True, project the gdf to latlong
-        latlong_crs = {'init':'epsg:4326'}
-        projected_gdf = gdf.to_crs(latlong_crs)
-        if not hasattr(gdf, 'gdf_name'):
-            gdf.gdf_name = 'unnamed'
-        log('Projected the GeoDataFrame "{}" to EPSG 4326 in {:,.2f} seconds'.format(gdf.gdf_name, time.time()-start_time))
-    else:
-        # else, project the gdf to UTM
-        # if GeoDataFrame is already in UTM, just return it
-        if (not gdf.crs is None) and ('proj' in gdf.crs) and (gdf.crs['proj'] == 'utm'):
-            return gdf
-
-        # calculate the centroid of the union of all the geometries in the GeoDataFrame
-        avg_longitude = gdf['geometry'].unary_union.centroid.x
+    # if gdf has no gdf_name attribute, create one now
+    if not hasattr(gdf, 'gdf_name'):
+        gdf.gdf_name = 'unnamed'
+
+    # if to_crs was passed-in, use this value to project the gdf
+    if to_crs is not None:
+        projected_gdf = gdf.to_crs(to_crs)
 
-        # calculate the UTM zone from this avg longitude and define the UTM CRS to project
-        utm_zone = int(math.floor((avg_longitude + 180) / 6.) + 1)
-        utm_crs = {'datum': 'NAD83',
-                   'ellps': 'GRS80',
-                   'proj' : 'utm',
-                   'zone' : utm_zone,
-                   'units': 'm'}
+    # if to_crs was not passed-in, calculate the centroid of the geometry to determine UTM zone
+    else:
+        if to_latlong:
+            # if to_latlong is True, project the gdf to latlong
+            latlong_crs = {'init':'epsg:4326'}
+            projected_gdf = gdf.to_crs(latlong_crs)
+            log('Projected the GeoDataFrame "{}" to EPSG 4326 in {:,.2f} seconds'.format(gdf.gdf_name, time.time()-start_time))
+        else:
+            # else, project the gdf to UTM
+            # if GeoDataFrame is already in UTM, just return it
+            if (not gdf.crs is None) and ('proj' in gdf.crs) and (gdf.crs['proj'] == 'utm'):
+                return gdf
+
+            # calculate the centroid of the union of all the geometries in the GeoDataFrame
+            avg_longitude = gdf['geometry'].unary_union.centroid.x
+
+            # calculate the UTM zone from this avg longitude and define the UTM CRS to project
+            utm_zone = int(math.floor((avg_longitude + 180) / 6.) + 1)
+            utm_crs = {'datum': 'NAD83',
+                       'ellps': 'GRS80',
+                       'proj' : 'utm',
+                       'zone' : utm_zone,
+                       'units': 'm'}
 
-        # project the GeoDataFrame to the UTM CRS
-        projected_gdf = gdf.to_crs(utm_crs)
-        if not hasattr(gdf, 'gdf_name'):
-            gdf.gdf_name = 'unnamed'
-        log('Projected the GeoDataFrame "{}" to UTM-{} in {:,.2f} seconds'.format(gdf.gdf_name, utm_zone, time.time()-start_time))
+            # project the GeoDataFrame to the UTM CRS
+            projected_gdf = gdf.to_crs(utm_crs)
+            log('Projected the GeoDataFrame "{}" to UTM-{} in {:,.2f} seconds'.format(gdf.gdf_name, utm_zone, time.time()-start_time))
 
     projected_gdf.gdf_name = gdf.gdf_name
     return projected_gdf
 
 
-def project_graph(G):
+def project_graph(G, to_crs=None):
     """
     Project a graph from lat-long to the UTM zone appropriate for its geographic location.
     
     Parameters
     ----------
     G : networkx multidigraph
-        the networkx graph to be projected to UTM
+        the networkx graph to be projected
+    to_crs : dict
+        if not None, just project to this CRS instead of to UTM
     
     Returns
     -------
@@ -128,7 +140,7 @@ def project_graph(G):
     log('Created a GeoDataFrame from graph in {:,.2f} seconds'.format(time.time()-start_time))
 
     # project the nodes GeoDataFrame to UTM
-    gdf_nodes_utm = project_gdf(gdf_nodes)
+    gdf_nodes_utm = project_gdf(gdf_nodes, to_crs=to_crs)
 
     # extract data for all edges that have geometry attribute
     edges_with_geom = []
@@ -142,7 +154,7 @@ def project_graph(G):
         gdf_edges = gpd.GeoDataFrame(edges_with_geom)
         gdf_edges.crs = G_proj.graph['crs']
         gdf_edges.gdf_name = '{}_edges'.format(G_proj.name)
-        gdf_edges_utm = project_gdf(gdf_edges)
+        gdf_edges_utm = project_gdf(gdf_edges, to_crs=to_crs)
 
     # extract projected x and y values from the nodes' geometry column
     start_time = time.time()

tests/test_osmnx.py

@@ -40,7 +40,7 @@ def test_imports():
 def test_gdf_shapefiles():
 
     city = ox.gdf_from_place('Manhattan, New York City, New York, USA')
-    city_projected = ox.project_gdf(city)
+    city_projected = ox.project_gdf(city, to_crs={'init':'epsg:3395'})
     ox.save_gdf_shapefile(city_projected)
 
     city = ox.gdf_from_place('Manhattan, New York City, New York, USA', buffer_dist=100)