基于Python与Geopandas围绕指定经纬度构建扇区网格的实现咨询

前置依赖 #

先安装所需第三方库：

pip install geopandas shapely pyproj numpy

实现逻辑 #

• 定义输入参数，包括中心点经纬度、扇区方位角、作用距离、张开角度、网格单元边长
• 调用大地测量工具计算扇区边界点位，生成闭合扇区面要素
• 基于扇区外接矩形生成规则网格，筛选出落在扇区范围内的网格单元
• 输出支持导出为shp、GeoJSON等常用地理数据格式

完整可运行代码 #

import geopandas as gpd
from shapely.geometry import Point, Polygon
from pyproj import Geod
import numpy as np

# ------------- 输入参数可自定义 -------------
center_lon = 116.4038  # 中心点经度
center_lat = 39.9151   # 中心点纬度
bearing = 90           # 扇区朝向方位角（0度为正北，顺时针递增）
distance = 5000        # 扇区作用距离，单位：米
sector_angle = 60      # 扇区张开角度，单位：度
grid_size = 500        # 单个网格单元边长，单位：米
# -------------------------------------------

# 初始化WGS84大地测量对象，用于方位、距离点位计算
geod = Geod(ellps='WGS84')
half_sector = sector_angle / 2

# 生成扇区面要素
left_az = bearing - half_sector
right_az = bearing + half_sector
# 采样圆弧点位保证扇区边界平滑
arc_points = []
for az in np.linspace(left_az, right_az, 100):
    lon, lat, _ = geod.fwd(center_lon, center_lat, az, distance)
    arc_points.append((lon, lat))
# 拼接扇区闭合边界
sector_coords = [(center_lon, center_lat)] + arc_points + [(center_lon, center_lat)]
sector_poly = Polygon(sector_coords)

# 坐标转投影坐标系保证网格长度精度，精度要求高可替换为对应区域的UTM投影
geo_crs = 'EPSG:4326'
proj_crs = 'EPSG:3857'
sector_gdf = gpd.GeoDataFrame(geometry=[sector_poly], crs=geo_crs).to_crs(proj_crs)
minx_proj, miny_proj, maxx_proj, maxy_proj = sector_gdf.total_bounds

# 生成规则网格
x_list = np.arange(minx_proj, maxx_proj + grid_size, grid_size)
y_list = np.arange(miny_proj, maxy_proj + grid_size, grid_size)
grid_list = []
for x in x_list[:-1]:
    for y in y_list[:-1]:
        grid = Polygon([(x, y), (x+grid_size, y), (x+grid_size, y+grid_size), (x, y+grid_size)])
        grid_list.append(grid)

# 筛选落在扇区内的网格，可将intersects替换为contains仅保留完全在扇区内的网格
grid_gdf = gpd.GeoDataFrame(geometry=grid_list, crs=proj_crs)
result_grid = grid_gdf[grid_gdf.intersects(sector_gdf.iloc[0].geometry)].to_crs(geo_crs)

# 结果导出（按需开启）
# result_grid.to_file('sector_grid.shp', encoding='utf-8')
# result_grid.to_file('sector_grid.geojson', driver='GeoJSON')

# 可视化验证（按需开启，需安装matplotlib）
# import matplotlib.pyplot as plt
# ax = sector_gdf.plot(alpha=0.2, color='red')
# result_grid.plot(ax=ax, edgecolor='gray', facecolor='none')
# plt.show()

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