OSMNX与Shapely切割多边形为网格仅得4块，代码问题排查求助

你调用ox.utils_geo._quadrat_cut_geometry时传入的quadrat_width参数单位与当前几何的坐标系单位不匹配：

• 你通过ox.geocode_to_gdf获取的巴黎市域几何默认使用WGS84坐标系（EPSG:4326），单位为十进制度
• 巴黎市域东西、南北跨度均仅为十余公里，对应十进制度单位下仅零点几度的差值
• 你设置的quadrat_width=750对应750度的切割宽度，远大于市域范围，因此仅能横竖各切一次，得到4个区块

先将几何投影到以米为单位的平面坐标系后再执行切割，示例修改代码如下：

import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import osmnx as ox
from descartes import PolygonPatch

# 获取巴黎市域几何
city = ox.geocode_to_gdf("Paris,France")
# 投影到法国官方平面坐标系EPSG:2154，单位为米
city_proj = ox.project_gdf(city, to_crs="EPSG:2154")
geometry = city_proj['geometry'].iloc[0]

# 此时quadrat_width=750对应750米宽的网格
geometry_split = ox.utils_geo._quadrat_cut_geometry(geometry, quadrat_width = 750)
polylist = [p for p in geometry_split]

# 绘图
west, south, east, north = city_proj.unary_union.bounds
fig, ax = plt.subplots(figsize=(30,30))
for polygon, n in zip(geometry_split, np.arange(len(polylist))):
    p = polygon.representative_point().coords[:][0]
    patch = PolygonPatch(polygon,ec='#100000',alpha=0.5, zorder=2)
    ax.add_patch(patch)
    plt.annotate(text=n, xy=p, horizontalalignment='center', size=15)
ax.set_xlim(west, east)
ax.set_ylim(south, north)
ax.axis('off')
plt.show()

# 如果需要将切割后的网格转回WGS84坐标系，可以执行以下操作
# split_gdf = gpd.GeoDataFrame(geometry=geometry_split, crs=city_proj.crs).to_crs(city.crs)

修改后即可得到按750米宽度切割的大量网格区块。

内容的提问来源于stack exchange，提问作者mark