Search code examples
geometrygeopandas

CRS the same but geometry in different numeric formats

I have two geometries with the same CRS (4326), but they have in are completely different X-Y axis formats. They should be overlapping. I am having issues with this geometry below. It says its 4326 but its the X/Y isn't in that projection.

import osmnx as ox
import networkx as nx
from shapely.geometry import Point, LineString, Polygon
import geopandas as gpd
from descartes import PolygonPatch

configure the place, network type, trip times, and travel speed

   place = 'Stockholm, Sweden'
    network_type = 'drive'
    trip_times = [15] #in minutes
    travel_speed = 4.5 #walking speed in km/hour
    G_4326 = ox.graph_from_place(place, network_type=network_type)

get nearest node

urban_intervention = ox.distance.get_nearest_node(G_4326, (59.33039855957031, 18.022981643676758))  
G_4326 = ox.project_graph(G_4326)

add an edge attribute for time in minutes required to traverse each edge

meters_per_minute = travel_speed * 1000 / 60 #km per hour to m per minute
for u, v, k, data in G_4326.edges(data=True, keys=True):
    data['time'] = data['length'] / meters_per_minute

isochrone_polys = []
for trip_time in sorted(trip_times, reverse=True):
    subgraph = nx.ego_graph(G_4326, urban_intervention, radius=trip_time, distance='time')
    node_points = [Point((data['x'], data['y'])) for node, data in subgraph.nodes(data=True)]
    bounding_poly = gpd.GeoSeries(node_points).unary_union.convex_hull
    isochrone_polys.append(bounding_poly)

convert to geopandas

treatment_radius = isochrone_polys[0]
treatment_radius_gdf = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[treatment_radius])
print(treatment_radius_gdf.crs)
print(treatment_radius_gdf)

result:

enter image description here

Any help would be very welcome!

Solution

  • You have a line of code that is projecting from EPSG:4326 to a UTM CRS. Additionally you are using very poor variable naming. you have called a variable G_4326 when it is not after projection!

    This is clearly documented: https://osmnx.readthedocs.io/en/stable/osmnx.html#osmnx.projection.project_graph

    # this line changes geometry from epsg:4326 to UTM CRS
G_4326 = ox.project_graph(G_4326)

    This makes sense given you are calculating in meters. However that means that the generated geometry CRS is not EPSG:4326. Set the CRS correctly, then you can project back to EPSG:4326. Clearly this means G_4326 is badly named, have not addressed that.

    output

    +proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +type=crs
                                            geometry
0  POLYGON ((676968.606 6569863.360, 676868.881 6...
epsg:4326
                                            geometry
0  POLYGON ((18.10176 59.23074, 18.10003 59.23086...

    full code

    import osmnx as ox
import networkx as nx
import geopandas as gpd
from shapely.geometry import Point

place = "Stockholm, Sweden"
network_type = "drive"
trip_times = [15]  # in minutes
travel_speed = 4.5  # walking speed in km/hour
G_4326 = ox.graph_from_place(place, network_type=network_type)
gdf_nodes1, gdf_edges1 = ox.graph_to_gdfs(G_4326)

# this line changes geometry from epsg:4326 to UTM CRS
G_4326 = ox.project_graph(G_4326)

gdf_nodes, gdf_edges = ox.graph_to_gdfs(G_4326)
m = gdf_edges.explore()

urban_intervention = ox.distance.get_nearest_node(
    G_4326, (59.33039855957031, 18.022981643676758)
)

meters_per_minute = travel_speed * 1000 / 60  # km per hour to m per minute
for u, v, k, data in G_4326.edges(data=True, keys=True):
    data["time"] = data["length"] / meters_per_minute

isochrone_polys = []
for trip_time in sorted(trip_times, reverse=True):
    subgraph = nx.ego_graph(
        G_4326, urban_intervention, radius=trip_time, distance="time"
    )
    node_points = [
        Point((data["x"], data["y"])) for node, data in subgraph.nodes(data=True)
    ]
    bounding_poly = gpd.GeoSeries(node_points).unary_union.convex_hull
    isochrone_polys.append(bounding_poly)
    
treatment_radius = isochrone_polys[0]
treatment_radius_gdf = gpd.GeoDataFrame(index=[0], crs=gdf_edges.crs, geometry=[treatment_radius])
print(treatment_radius_gdf.crs)
print(treatment_radius_gdf)
print(treatment_radius_gdf.to_crs("epsg:4326").crs)
print(treatment_radius_gdf.to_crs("epsg:4326"))
treatment_radius_gdf.explore(m=m, color="red")