一种用 Python 创建鱼网格（fishnet grid）的方法！

这篇文章介绍了一种用 Python 创建鱼网格（fishnet grid）的方法：先在目标多边形内部生成等间距点阵，再围绕这些点生成 Voronoi/Thiessen 多边形，最后把这些多边形裁剪回原始边界。这样得到的网格在内部区域可以保持单元面积相等，边界处则会因为被原始区域裁剪而变成不完整单元。

作者的场景是地理空间分析：把连续的城市空间离散化成规则网格，便于按网格聚合坐标点、计算统计指标，并进一步叠加其他空间数据。

如果你想在一个多边形区域内创建单元大小几乎完全一致的鱼网格（边界处被裁剪的格子除外），可以先在多边形内部生成一组等间距点，再围绕每个点生成 Thiessen 多边形，也就是 Voronoi 多边形。本文末尾给出了完整代码。所有网格生而平等，但有些网格比其他网格更平等。

为什么要做网格

在做地理空间数据分析时，最常见的第一步之一，就是把多个坐标点聚合到研究区域上的一个统一网格里。这样做的好处很直接：你可以围绕采集到的数据计算各种统计指标，也可以继续往网格里叠加更多数据，让模型获得更丰富的空间特征。从数据处理的角度看，这其实就是地理空间版的“把连续变量离散化”或者“分箱”。如果你随手 Google 一下，或者问问 ChatGPT，会看到许多许多创建网格或鱼网格的方法。本文不是一篇比较各种方法优劣的论文，而只是记录我在项目中为了做出合适网格时，偶然摸索到的一种办法。下面开始。

选择并可视化网格区域

首先，导入这个项目需要用到的全部 Python 包。

import numpy as npimport geopandas as gpdimport foliumimport warnings from shapely.prepared import prepfrom shapely.geometry import Point,Polygonfrom shapely.ops import cascaded_unionfrom scipy import spatial warnings.filterwarnings('ignore')

以某市里的几个区作为测试案例。
读取 shapefile，并筛选出需要的行政区file_path = '../data/district_shapefiles/District_Boundary.shp'districts_gpd = gpd.read_file(file_path)khi_districts_list = ['KARACHI SOUTH','KARACHI EAST','CLIFTON CANTONMENT']khi_districts_gpd = districts_gpd[districts_gpd.DISTRICT</span>.isin(khi_districts_list)]khi_districts_gpd = khi_districts_gpd.reset_index(drop=True)

接下来，我们把三个区域合并成一个多边形。这个合并后的多边形，就是后面整个网格需要覆盖的总区域。这里用 Folium 把它可视化出来。
使用 cascaded union 将多个多边形合并为一个khi_shape = cascaded_union(khi_districts_gpd.geometry) loc = list(khi_districts_gpd.geometry[0].exterior.coords)[0]loc = loc[::-1]map_ =  folium.Map(location=loc,zoom_start=12)map_.add_child(folium.GeoJson(data=khi_shape))map_

在选定多边形内生成点阵

既然我们已经拿到了表示总覆盖区域的 polygon 对象，下一步就是在这个多边形内部生成一组等距点。这里有一个细节很重要：经纬度坐标的单位是角度，不适合直接拿来计算“米”级距离。因此代码会先把坐标系从 WGS84 转成 Web Mercator，这样点与点之间的距离就可以按米来计算，而不是按弧度或角度来算。
将投影从 WGS84 转换到 Web Mercator# 这样点阵之间的距离就可以按“米”来计算，而不是按经纬度角度计算 poly = gpd.GeoDataFrame({'geometry':[khi_shape]},crs=4326)</span>.to_crs(epsg=3857)['geometry'].values[0] # 获取选定多边形的外接矩形边界latmin, lonmin, latmax, lonmax = poly.boundsprep_polygon = prep(poly)valid_points = []points = [] # 点阵中相邻两个点之间的距离为 500 米# 你可以调整这个值来控制网格单元大小resolution = 500 # 创建点阵for lat in np.arange(latmin, latmax, resolution):    for lon in np.arange(lonmin, lonmax, resolution):        points.append(Point((round(lat,4), round(lon,4)))) # 只保留落在选定多边形内部的点valid_points.extend(filter(prep_polygon.contains,points))# 将投影转回 WGS84，并把点列表转换为 GeoDataFramemesh_gpd = gpd.GeoDataFrame({'geometry':valid_points},crs=3857)</span>.to_crs(epsg=4326) mesh_gpd = mesh_gpd.reset_index()我们再用 Folium 把这个点阵可视化一下。map_ =  folium.Map(location=loc,zoom_start=12)for ind, row in mesh_gpd.iterrows():    point = row.geometry    lat, lon = point.y, point.x    folium.CircleMarker([lat, lon], radius=0.1,                          color='blue', fill=True,                          fill_color='blue').add_to(map_) map_

生成 Voronoi 剖分

这些点阵中的点，会作为 Voronoi 剖分的中心点。不过继续写代码之前，我们先稍微绕开一下，解释一下什么是 Voronoi Tessellation，也就是 Thiessen Polygon。用数学语言来说，Voronoi 图是把一个平面按照一组给定对象划分成若干区域。最简单的情况是：这些对象就是平面上的有限多个点，它们通常被称为 seed、site 或 generator。对于每个种子点，都会有一个对应区域，叫 Voronoi cell。这个区域包含了平面上所有“离这个种子点比离其他种子点都更近”的点。换句话说，如果一个随机点出现在某个 seed 或 site 的区域内，那么这个区域的 seed 就是离它最近的那个点。Towards Data Science 上有一篇很棒的文章更深入地介绍了这个概念，如果你想继续了解，可以去读一读。原文中的示意图也来自那里。为了创建我们的网格，需要围绕前面生成的等距点阵创建 Voronoi 多边形。因为这些种子点之间的距离是相等的，所以在选定多边形范围内，最终得到的区域或多边形大小也会一致。现在运行最后一段代码，生成我们的网格。
从 geometry 对象中提取 x、y 坐标，并用它们生成 Voronoi 剖分mesh_gpd['lng_lat'] = mesh_gpd.geometry</span>.apply(lambda coord: (round(coord.x,5),round(coord.y,5)))vor = spatial.Voronoi(list(mesh_gpd['lng_lat'].values))# 这个函数来自 GitHub 上某个 gist，用于把边界处无限延伸的区域转换成有限区域。# 它能正常工作。抱歉，我现在已经找不到原始 gist，# 所以无法给原作者补上应有的署名。def func_voronoi_finite_polygons_2d(vor, radius=None):    """    将二维Voronoi 图中的无限区域重构为有限区域。    参数    ----------    vor: Voronoi        输入的Voronoi 图。    radius: float, optional        “无穷远点”的距离。    返回    -------    regions: list of tuples        每个修正后Voronoi 区域中顶点的索引。    vertices: list of tuples        修正后Voronoi 顶点的坐标。它等于输入顶点坐标，        并在末尾追加了“无穷远点”。    """    ifvor.points.shape[1] != 2:        raiseValueError("Requires 2D input")    new_regions= []    new_vertices= vor.vertices.tolist()    center= vor.points.mean(axis=0)    ifradius is None:        radius= vor.points.ptp().max()*2    #为每个点构造一张包含所有 ridge 的映射    all_ridges= {}    for(p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):        all_ridges.setdefault(p1,[]).append((p2, v1, v2))        all_ridges.setdefault(p2,[]).append((p1, v1, v2))    #重构无限区域    forp1, region in enumerate(vor.point_region):        vertices= vor.regions[region]        ifall(v >= 0 for v in vertices):            #有限区域            new_regions.append(vertices)            continue        #重构非有限区域        ridges= all_ridges[p1]        new_region= [v for v in vertices if v >= 0]        forp2, v1, v2 in ridges:            ifv2 < 0:                v1,v2 = v2, v1            ifv1 >= 0:                #有限 ridge：已经在区域中                continue            #计算无限 ridge 缺失的端点            t= vor.points[p2] - vor.points[p1] # 切向量            t/= np.linalg.norm(t)            n= np.array([-t[1], t[0]])# 法向量            midpoint= vor.points[[p1, p2]].mean(axis=0)            direction= np.sign(np.dot(midpoint - center, n)) * n            far_point= vor.vertices[v2] + direction * radius            new_region.append(len(new_vertices))            new_vertices.append(far_point.tolist())        #将区域顶点按逆时针排序        vs= np.asarray([new_vertices[v] for v in new_region])        c= vs.mean(axis=0)        angles= np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])        new_region= np.array(new_region)[np.argsort(angles)]        #完成        new_regions.append(new_region.tolist())    returnnew_regions, np.asarray(new_vertices)# 将边界处的无限区域转换为有限区域regions,vertices=func_voronoi_finite_polygons_2d(vor,radius=0.07)# 这个函数用于把结果多边形裁剪到原始多边形边界内def func_intersect_vor(row):    returnkhi_shape.intersection(row.vor)# 将 regions 和 vertices 转换为 Polygon 对象，# 并得到最终网格对应的 GeoDataFramepolys = []for region in regions:    poly= vertices[region]    polys.append(Polygon(poly))vor_gpd = gpd.GeoDataFrame(polys)vor_gpd = vor_gpd.rename(columns={    0:'vor'})vor_gpd['bounded_vor'] = vor_gpd.apply(func_intersect_vor,1)vor_gpd = gpd.GeoDataFrame(vor_gpd,geometry='bounded_vor')vor_gpd = vor_gpd.drop(columns='vor')vor_gpd = vor_gpd.reset_index()vor_gpd = vor_gpd.rename(columns={'index':'vor_id'})vor_gpd['vor_id'] = vor_gpd.vor_id.apply(lambda x: 'vor'+str(x))我们的网格已经准备好了。再用 Folium 把它可视化出来。map_ =folium.Map(location=loc,zoom_start=12)for ind, row in vor_gpd.iterrows():    geometry_obj= row.bounded_vor    map_.add_child(folium.GeoJson(geometry_obj)) map_

最终结果

这里的 vor_gpd 变量是一个 GeoDataFrame。它的 geometry 列里，每一个 polygon 对象都是鱼网格中的一个单元。换个角度看，它也可以被理解为对应 Voronoi 剖分中某个种子点的区域。

这套方法什么时候适合用

如果你的目标是把城市空间、行政区域、服务覆盖范围或任意不规则区域切成相对均匀的小单元，这种“等距点阵 + Voronoi + 边界裁剪”的思路非常实用。它比简单画矩形格子更贴合边界，也比手工分区更容易自动化。不过，如果你要做严肃的面积统计，建议根据研究区域选择更合适的本地投影，而不是直接依赖 Web Mercator。边界处被裁剪的小格子也最好单独标记，后续统计时可以按面积做修正。

以上就是“一种用 Python 创建鱼网格（fishnet grid）的方法！”的详细内容，想要了解更多Python教程欢迎持续关注编程学习网。