当前位置：首页>python>一种用 Python 创建鱼网格(fishnet grid)的方法

一种用 Python 创建鱼网格(fishnet grid)的方法

2026-06-18 06:06:57

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

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

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

所有网格生而平等，但有些网格比其他网格更平等。

为什么要做网格

在做地理空间数据分析时，最常见的第一步之一，就是把多个坐标点聚合到研究区域上的一个统一网格里。这样做的好处很直接：你可以围绕采集到的数据计算各种统计指标，也可以继续往网格里叠加更多数据，让模型获得更丰富的空间特征。

从数据处理的角度看，这其实就是地理空间版的“把连续变量离散化”或者“分箱”。

如果你随手 Google 一下，或者问问 ChatGPT，会看到许多许多创建网格或鱼网格的方法。本文不是一篇比较各种方法优劣的论文，而只是记录我在项目中为了做出合适网格时，偶然摸索到的一种办法。下面开始。

选择并可视化网格区域

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

import numpy as npimport geopandas as gpdimport foliumimport warningsfrom shapely.prepared import prepfrom shapely.geometry import Point,Polygonfrom shapely.ops import cascaded_unionfrom scipy import spatialwarnings.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):forlon 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.geometrylat,lon = point.y, point.xfolium.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, v1ifv1 >= 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)) * nfar_point= vor.vertices[v2] + direction * radiusnew_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_vormap_.add_child(folium.GeoJson(geometry_obj))map_