如何用geopandas实现GIS地图网格划分？

Question

如何使用GeoPandas基于给定地理范围生成等间距的矢量网格（如1km×1km的正方形网格）并保存为GIS格式？常见问题包括：如何从Shapefile或GeoJSON读取边界，利用边界最小外接矩形（bounds）计算网格坐标，结合Shapely构造多边形网格，再通过GeoDataFrame组织数据。难点在于坐标参考系统（CRS）的选择——需先投影到合适的投影坐标系（如UTM），以确保网格距离准确，避免在经纬度坐标系下因变形导致网格不均匀。最后，如何将生成的网格导出为Shapefile或GeoJSON以便在QGIS等软件中可视化？

Answer 1

一、引言：地理网格在空间分析中的核心作用

在GIS与空间数据分析中，生成规则的矢量网格（如1km×1km正方形）是进行区域划分、热力图可视化、空间聚合统计和模型输入的基础操作。GeoPandas作为Python生态中处理矢量地理数据的核心库，结合Shapely、Fiona等工具，能够高效实现此类任务。然而，实际应用中常因坐标参考系统（CRS）选择不当导致网格变形，影响精度。本文将从基础到进阶，系统阐述如何使用GeoPandas生成等间距矢量网格，并导出为标准GIS格式。

二、数据准备与边界读取

首先需加载指定区域的边界数据，支持Shapefile或GeoJSON格式。GeoPandas提供统一接口gpd.read_file()进行读取：

import geopandas as gpd

# 读取边界文件
boundary = gpd.read_file("data/boundary.shp")  # 或 "boundary.geojson"
print(boundary.crs)  # 查看原始CRS
    

输出结果通常包含几何列（geometry）及属性字段，其CRS信息至关重要，直接影响后续投影与距离计算准确性。

三、坐标参考系统（CRS）转换的重要性

地理坐标系（如WGS84, EPSG:4326）以经纬度表示，单位为度，不适合距离测量。直接在此基础上生成1km网格会导致高纬度地区严重压缩。因此必须先投影至合适的投影坐标系，推荐使用UTM分区（Universal Transverse Mercator），其单位为米，局部形变小。

自动确定最佳UTM带的函数如下：

def get_utm_crs(gdf):
    """根据GeoDataFrame中心点自动获取UTM CRS"""
    centroid = gdf.geometry.unary_union.centroid
    lon, lat = centroid.x, centroid.y
    zone = int((lon + 180) // 6) + 1
    epsg = 32600 + zone if lat >= 0 else 32700 + zone
    return f"EPSG:{epsg}"

# 投影到UTM
utm_crs = get_utm_crs(boundary)
boundary_proj = boundary.to_crs(utm_crs)
    

四、基于最小外接矩形生成网格坐标

利用投影后的边界获取最小外接矩形（bounds），并按设定分辨率（如1000米）生成网格点阵：

from shapely.geometry import box
import numpy as np

def create_grid(gdf, cell_size=1000):
    bounds = gdf.total_bounds  # [minx, miny, maxx, maxy]
    minx, miny, maxx, maxy = bounds
    
    # 向下/左扩展至最近的cell_size倍数
    x_steps = np.arange(np.floor(minx / cell_size) * cell_size, 
                        np.ceil(maxx / cell_size) * cell_size, 
                        cell_size)
    y_steps = np.arange(np.floor(miny / cell_size) * cell_size, 
                        np.ceil(maxy / cell_size) * cell_size, 
                        cell_size)

    grid_polygons = []
    for i in range(len(x_steps)-1):
        for j in range(len(y_steps)-1):
            poly = box(x_steps[i], y_steps[j], 
                      x_steps[i+1], y_steps[j+1])
            grid_polygons.append(poly)
    
    return grid_polygons

grid_cells = create_grid(boundary_proj, cell_size=1000)
    

五、构建GeoDataFrame并保留原始CRS信息

将生成的多边形列表封装为GeoDataFrame，并恢复原始地理坐标系以便与其他数据集成：

grid_gdf = gpd.GeoDataFrame({'geometry': grid_cells}, crs=utm_crs)
# 可选：裁剪网格至边界范围内
grid_clipped = gpd.overlay(grid_gdf, boundary_proj, how='intersection')

# 转回地理坐标系（便于发布）
grid_final = grid_clipped.to_crs(boundary.crs)
    

六、导出为标准GIS格式供QGIS等软件使用

支持多种格式导出，常用包括Shapefile和GeoJSON：

格式	优点	缺点	代码示例
Shapefile	兼容性强，QGIS/ArcGIS通用	不支持嵌套属性，编码问题常见	grid_final.to_file("output/grid.shp")
GeoJSON	文本可读，Web友好	大文件性能差	grid_final.to_file("output/grid.geojson", driver="GeoJSON")
GeoPackage	单文件多图层，现代标准	旧软件可能不支持	grid_final.to_file("output/grid.gpkg", driver="GPKG")

七、完整流程的Mermaid流程图

graph TD A[读取边界文件
(Shapefile/GeoJSON)] --> B{检查CRS} B -- 地理坐标系 --> C[投影至UTM] B -- 投影坐标系 --> D[继续处理] C --> D D --> E[计算最小外接矩形] E --> F[生成等距网格坐标] F --> G[构造Shapely多边形] G --> H[创建GeoDataFrame] H --> I[可选：裁剪至边界内] I --> J[转回原始CRS] J --> K[导出为Shapefile/GeoJSON] K --> L[在QGIS中可视化验证]

八、常见问题与调试建议

• 网格间距不准？ 检查是否未投影至米制坐标系。
• 导出后属性丢失？ Shapefile对字段名长度和类型敏感，建议使用GeoPackage。
• 性能瓶颈？ 大范围高密度网格可用Dask+GeoPandas分块处理。
• 跨UTM带怎么办？ 对于超大区域，应分带处理或采用Albers等区域投影。
• 如何添加唯一ID？ 在GeoDataFrame中增加索引列：grid_final['grid_id'] = range(len(grid_final))

九、扩展应用场景

该方法不仅限于正方形网格，还可拓展为：

1. 六边形蜂窝网格（使用h3-py或自定义算法）
2. 动态分辨率网格（基于人口密度调整cell_size）
3. 三维体素网格（结合DEM高程数据）
4. 时空立方体构建（叠加时间维度）
5. 作为机器学习的空间特征容器（每个网格存储统计特征）
6. 用于无人机航拍路径规划的覆盖网格
7. 应急响应中的责任区划分
8. 城市功能区识别的基础单元
9. 生态监测样方布设辅助工具
10. 电信基站信号覆盖模拟网格