CRS 相同但几何形状不同的数字格式

Question

我有两个具有相同 CRS (4326) 的几何图形，但它们的 X-Y 轴格式完全不同。它们应该重叠。我对下面的这个几何有问题。它说它是 4326，但它的 X/Y 不在那个投影中。

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

配置地点、网络类型、行程时间和行驶速度

   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)

获取最近的节点

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

为遍历每条边所需的时间添加边属性（以分钟为单位）

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)

转换为 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)

结果：

非常欢迎任何帮助！

Answer 1

您有一行代码从 EPSG:4326 投影到 UTM CRS。此外，您使用的变量命名非常糟糕。你调用了一个变量 G_4326 而不是在投影之后！

这已明确记录：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)

考虑到您以米为单位计算，这是有道理的。然而，这意味着生成的几何 CRS 不是 EPSG:4326。正确设置CRS，则可以投影回EPSG:4326。显然，这意味着 G_4326 命名错误，尚未解决。

输出

+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...

完整代码

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")

CRS 相同但几何形状不同的数字格式

CRS the same but geometry in different numeric formats

geometry

geopandas

输出

完整代码