渔网网格，然后在每个网格点周围创建泰森多边形。完整的代码可以在本文的末尾找到。
所有网格都是平等的，有些网格比其他网格更平等。
分析地理空间数据的第一步是以某种方式将多个坐标数据点聚合在要分析的区域上的统一网格上。这将帮助您计算与您正在收集的数据相关的不同统计数据，同时还允许您通过向网格添加额外的数据点来进一步丰富您的模型。这与离散连续变量或分箱数据的地理空间等效。
如果您 Google 搜索（或 chatGTP 上的提示），您将看到许多创建网格或鱼网网格的方法。这篇文章并不是关于哪种技术更好或更坏的论文，而是试图记录我在尝试为我的项目制作正确的网格时偶然发现的一种技术。所以就这样了。
选择并可视化网格区域。
首先，我们导入该项目所需的所有包。
mport numpy as np
import geopandas as gpd
import folium
import warnings
from shapely.prepared import prep
from shapely.geometry import Point,Polygon
from shapely.ops import cascaded_union
from scipy import spatial
warnings.filterwarnings('ignore')
# read and filter relevant shapefiles 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\ .isin(khi_districts_list)] khi_districts_gpd = khi_districts_gpd.reset_index(drop=True)
接下来，我们使用 folium 来可视化这一点
# Combine multiple polygons into one using a 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_
在选定的多边形内生成点网格
由于我们拥有多边形对象中的总覆盖区域，因此下一步是在该多边形内创建等距点的网格。
# Changing projection form WGS84 to Web Mercator for caclulating distance
# between mesh points in meters not radians.
poly = gpd.GeoDataFrame({'geometry':[khi_shape]},crs=4326)\
.to_crs(epsg=3857)['geometry'].values[0]
#Get the bounding points of the selected polygon.
latmin, lonmin, latmax, lonmax = poly.bounds
prep_polygon = prep(poly)
valid_points = []
points = []
# Each point in mesh will be 500 meters apart.
# You can play with this to adjust the cell size
resolution = 500
# create mesh
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))))
# only keep those points that are within the the selected polygon.
valid_points.extend(filter(prep_polygon.contains,points))
# change projection back to WGS84 convert list of points to Geodataframe.
mesh_gpd = gpd.GeoDataFrame({'geometry':valid_points},crs=3857)\
.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 曲面细分或泰森多边形。
用数学术语来说，Voronoi 图是将平面划分为靠近给定对象组中每个对象的区域。在最简单的情况下，这些对象只是平面上有限的多个点（称为种子、站点或生成器）。对于每个种子，都有一个相应的区域，称为 Voronoi 单元，由平面上距离该种子比距离其他任何位置更近的所有点组成。
本质上，如果随机点在现有种子或站点的区域内生成，则它是最接近该区域的种子。关于数据科学有一篇很棒的博客，它更详细地探讨了这个概念，如果您想了解更多信息，我建议您阅读它。下面的可视化是从同一个地方偷来的。
为了创建网格，我们需要围绕我们创建的等距网格创建 Voronoi 镶嵌。由于种子的间距相等，因此生成的区域或多边形在我们所选多边形的边界内具有相同的大小。
现在让我们运行最后一段代码来获取网格
# Getting x and y coordinates from geometry object and using it to
# generate tesselations.
mesh_gpd['lng_lat'] = mesh_gpd.geometry\
.apply(lambda coord: (round(coord.x,5),round(coord.y,5)))
vor = spatial.Voronoi(list(mesh_gpd['lng_lat'].values))
# Function I got off some gist from GitHub that converts infinite
# bordering regions to finite regions.It just works? Apologies to
# the original author, I am unable to locate the original gist and
# give due credit.
def func_voronoi_finite_polygons_2d(vor, radius=None):
"""
Reconstruct infinite voronoi regions in a 2D diagram to finite
regions.
Parameters
----------
vor : Voronoi
Input diagram
radius : float, optional
Distance to 'points at infinity'.
Returns
-------
regions : list of tuples
Indices of vertices in each revised Voronoi regions.
vertices : list of tuples
Coordinates for revised Voronoi vertices. Same as coordinates
of input vertices, with 'points at infinity' appended to the
end.
"""
if vor.points.shape[1] != 2:
raise ValueError("Requires 2D input")
new_regions = []
new_vertices = vor.vertices.tolist()
center = vor.points.mean(axis=0)
if radius is None:
radius = vor.points.ptp().max()*2
# Construct a map containing all ridges for a given point
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))
# Reconstruct infinite regions
for p1, region in enumerate(vor.point_region):
vertices = vor.regions[region]
if all(v >= 0 for v in vertices):
# finite region
new_regions.append(vertices)
continue
# reconstruct a non-finite region
ridges = all_ridges[p1]
new_region = [v for v in vertices if v >= 0]
for p2, v1, v2 in ridges:
if v2 < 0:
v1, v2 = v2, v1
if v1 >= 0:
# finite ridge: already in the region
continue
# Compute the missing endpoint of an infinite ridge
t = vor.points[p2] - vor.points[p1] # tangent
t /= np.linalg.norm(t)
n = np.array([-t[1], t[0]]) # normal
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())
# sort region counterclockwise
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)]
# finish
new_regions.append(new_region.tolist())
return new_regions, np.asarray(new_vertices)
# Convering infinite border regions to finite regions.
regions,vertices=func_voronoi_finite_polygons_2d(vor,radius=0.07)
# Funciton that clips our resulting polygon within the
# bounds of our original polygon.
def func_intersect_vor(row):
return khi_shape.intersection(row.vor)
# Converting regions and vertices to Polygon objects
# and getting a geodataframe with the final grid.
polys = []
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，其中几何列中的每个多边形对象都是渔网网格的一个单元格，或者是 Voroni 曲面细分的相应种子的一个区域 - 取决于您想要如何查看它。