geopandas：生成研究区正方形格网、正六边形格网和同心圆环

本文章均以经纬度为单位进行分割
如果希望用m为单位，给地图加个投影就行

一、正方形格网

def make_grid(boundary, radius):
	import geopandas as gpd
	import numpy as np
	from shapely.geometry import Point, Polygon
	
    # 以边界的中心点按照radius半径进行外扩
    minx, miny, maxx, maxy = boundary.total_bounds
    deltax = ((maxx-minx)%(radius*2)) / 2
    deltay = ((maxy-miny)%(radius*2)) / 2
    minx -= deltax
    maxx += deltax
    miny -= deltay
    maxy += deltay
    cells = []
    for x in np.arange(minx, maxx, radius*2):
        for y in np.arange(miny, maxy, radius*2):
            cells.append(shapely.geometry.Polygon([(x,y),(x,y+radius*2),(x+radius*2,y+radius*2),(x+radius*2,y),(x,y)]))
    return gpd.GeoDataFrame(cells,columns=['geometry'],crs=boundary.crs).sjoin(boundary,how='inner')  

import matplotlib.pyplot as plt

city = gpd.read_file('your data', encoding='gbk')
gdf = make_grid(city, 0.05)

fig, ax = plt.subplots(1,2,figsize=(20,10))
ax[0] = city.plot(ax=ax[0])
ax[1] = gdf.boundary.plot(ax=ax[1])

二、正六边形格网

def make_hexagon(gdf, radius):
	import geopandas as gpd
	import numpy as np
	import math
	from shapely.geometry import Point, Polygon
	
    # 以边界的中心点按照radius半径进行外扩
    minx, miny, maxx, maxy = gdf.total_bounds
    deltax = ((maxx-minx)%(radius*2))/2
    deltay = ((maxy-miny)%(radius*math.sqrt(3)))/2
    minx -= deltax
    maxx += deltax
    miny -= deltay
    maxy += deltay
    cells = []
    for x in np.arange(minx, maxx, radius*3):
        for y in np.arange(miny, maxy, radius*math.sqrt(3)):
            cells.append(Polygon([(x,y),(x-radius/2,y+radius*math.sqrt(3)/2),(x,y+radius*math.sqrt(3)),(x+radius,y+radius*math.sqrt(3)),(x+1.5*radius,y+radius*math.sqrt(3)/2),(x+radius,y),(x,y)]))
    for x in np.arange(minx+1.5*radius,maxx,radius*3):
        for y in np.arange(miny+radius*math.sqrt(3)/2,maxy,radius*math.sqrt(3)):
            cells.append(Polygon([(x,y),(x-radius/2,y+radius*math.sqrt(3)/2),(x,y+radius*math.sqrt(3)),(x+radius,y+radius*math.sqrt(3)),(x+1.5*radius,y+radius*math.sqrt(3)/2),(x+radius,y),(x,y)]))
    return gpd.GeoDataFrame(cells,columns=['geometry'],crs=gdf.crs).sjoin(gdf,how='inner')

import matplotlib.pyplot as plt

city = gpd.read_file('your data', encoding='gbk')
gdf = make_hexagon(city, 0.05)

fig, ax = plt.subplots(1,2,figsize=(20,10))
ax[0] = city.plot(ax=ax[0])
ax[1] = gdf.boundary.plot(ax=ax[1])

三、同心环

def make_rings(gdf, widths):
    # 根据凸包确定最长半径
    boundary = gdf.unary_union
    convex_hulls = boundary.convex_hull  # 计算凸包
    convex_hull_coords = convex_hulls.exterior.coords  # 获取凸包的顶点坐标
    max_distance = 0
    for i in range(len(convex_hull_coords)):  # 计算每对顶点之间的距离，并找到最长距离
        for j in range(len(convex_hull_coords)):
            point1 = Point(convex_hull_coords[i])
            point2 = Point(convex_hull_coords[j])
            distance = point1.distance(point2)
            if distance > max_distance:
                max_distance = distance             
    # 获取城市中心的经度和纬度
    citycenter = Point(convex_hulls.centroid)  
    # 根据凸包确定的半径寻找包裹全部点的半径
    R = max_distance / 2
    while True:
        all_points_inside = True
        for point in convex_hull_coords:
            point = Point(point)
            distance_to_center = point.distance(citycenter)
            if distance_to_center > R:
                all_points_inside = False
                break
        if all_points_inside:
            break
        else:
            R += widths  # 增加1个单位
    # 最终的环数目
    # num_rings = int(R // widths) + 1
    num_rings = int(R // widths)
    # 创建同心环的几何对象
    rings = []
    for i in range(num_rings):  # 创建指定数量的同心环
        if i == 0:
            radius = widths
            ring = citycenter.buffer(radius)
        else:
            inner_radius = i * widths  # 内半径
            outer_radius = (i + 1) * widths  # 外半径
            # 生成环
            ring = citycenter.buffer(outer_radius) - citycenter.buffer(inner_radius)
        rings.append(ring)
    return gpd.GeoDataFrame(data={'ID': range(len(rings)), 'Width': range(1, num_rings + 1)}, geometry=rings, crs=gdf.crs)

city = gpd.read_file('your data', encoding='gbk')
gdf = make_rings(city, 0.1)

ax = gdf.boundary.plot(figsize=(10,10),color='k')
city.plot(ax=ax)