Python应用指南：获取行政区最小外接矩形

有的时候我们用行政区的面层进行裁剪的时候，特别是在基于行政区裁剪的渔网，因为行政边界本身就是不规则的原因，在边界的渔网匹配数据会出现匹配不上的问题，其实如果周边有数据完全可以采用最近匹配的原则，但是就是因为行政区裁太可丁可卯了，会丢失一部分数据，这个时候可以通过形成一个最小外接矩形的形式，可以让数据裁剪保持一个微妙的范围，不至于太大，又不至于没有二次处理的空间，第二个方面就是比如我们需要行政区里面的POI，以高德为例，因为行政区本身不是规则的，我们通过建立外接矩阵，这样获取其中的数据就不会出现丢失的情况了；

行政区边界获取部分，这里采用了阿里云旗下的数据可视化平台，其数据源本身也是来着高德地图，数据质量还是比较可靠的；

全国行政区划边界GeoJSON数据下载平台阿里云：DataV.GeoAtlas地理小工具系列 (aliyun.com)

这里把下载到本地的geojson或者json的文件路径改成自己的即可；

完整代码#运行环境Python 3.11

python">import geopandas as gp# 结果变量
resultMinLon = 99999.0  # 最小经度
resultMaxLon = -99999.0  # 最大经度
resultMinLat = 99999.0  # 最小纬度
resultMaxLat = -99999.0  # 最大纬度# 读取厦门市边界数据
# 注意：请确保文件路径正确
loadGeoData = gp.read_file("D:\\data\\厦门市.json")# 遍历每一行数据
for i in range(len(loadGeoData)):# 读取几何数据loadGeometry = loadGeoData.loc[i, "geometry"]# 检查几何数据类型，可能是单个多边形（Polygon）或多个多边形（MultiPolygon）if loadGeometry.geom_type == 'Polygon':# 单个多边形# 遍历多边形的所有顶点for z in range(len(loadGeometry.exterior.coords)):# 获取当前顶点的经度和纬度loadLon, loadLat = loadGeometry.exterior.coords[z]# 更新最小经度if loadLon < resultMinLon:resultMinLon = loadLon# 更新最大经度if loadLon > resultMaxLon:resultMaxLon = loadLon# 更新最小纬度if loadLat < resultMinLat:resultMinLat = loadLat# 更新最大纬度if loadLat > resultMaxLat:resultMaxLat = loadLatelif loadGeometry.geom_type == 'MultiPolygon':# 多个多边形# 遍历每个多边形for polygon in loadGeometry.geoms:# 遍历多边形的所有顶点for z in range(len(polygon.exterior.coords)):# 获取当前顶点的经度和纬度loadLon, loadLat = polygon.exterior.coords[z]# 更新最小经度if loadLon < resultMinLon:resultMinLon = loadLon# 更新最大经度if loadLon > resultMaxLon:resultMaxLon = loadLon# 更新最小纬度if loadLat < resultMinLat:resultMinLat = loadLat# 更新最大纬度if loadLat > resultMaxLat:resultMaxLat = loadLat# 结果输出
print("厦门市最小经度：", resultMinLon)
print("厦门市最大经度：", resultMaxLon)
print("厦门市最小纬度：", resultMinLat)
print("厦门市最大纬度：", resultMaxLat)

这里会直接打印出来两个对角线坐标，另外，因为数据源是高德的，用的是（GCJ-02）坐标系，最好先转成WGS84地理坐标系或者CGCS2000地理坐标系。如果需要可视化的话,可以跑一下下面的代码，与上面的区别就是通过调用了matplotlib包，使得结果可视化；

完整代码#运行环境Python 3.11

python">import geopandas as gp
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle# 配置 matplotlib 支持中文显示
plt.rcParams['font.family'] = ['SimHei']  # 使用黑体字体
plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号# 结果变量
resultMinLon = 99999.0  # 最小经度
resultMaxLon = -99999.0  # 最大经度
resultMinLat = 99999.0  # 最小纬度
resultMaxLat = -99999.0  # 最大纬度# 读取厦门市边界数据
# 注意：请确保文件路径正确
loadGeoData = gp.read_file("D:\\data\\厦门市.json")# 遍历每一行数据
for i in range(len(loadGeoData)):# 读取几何数据loadGeometry = loadGeoData.loc[i, "geometry"]# 检查几何数据类型，可能是单个多边形（Polygon）或多个多边形（MultiPolygon）if loadGeometry.geom_type == 'Polygon':# 单个多边形# 遍历多边形的所有顶点for z in range(len(loadGeometry.exterior.coords)):# 获取当前顶点的经度和纬度loadLon, loadLat = loadGeometry.exterior.coords[z]# 更新最小经度if loadLon < resultMinLon:resultMinLon = loadLon# 更新最大经度if loadLon > resultMaxLon:resultMaxLon = loadLon# 更新最小纬度if loadLat < resultMinLat:resultMinLat = loadLat# 更新最大纬度if loadLat > resultMaxLat:resultMaxLat = loadLatelif loadGeometry.geom_type == 'MultiPolygon':# 多个多边形# 遍历每个多边形for polygon in loadGeometry.geoms:# 遍历多边形的所有顶点for z in range(len(polygon.exterior.coords)):# 获取当前顶点的经度和纬度loadLon, loadLat = polygon.exterior.coords[z]# 更新最小经度if loadLon < resultMinLon:resultMinLon = loadLon# 更新最大经度if loadLon > resultMaxLon:resultMaxLon = loadLon# 更新最小纬度if loadLat < resultMinLat:resultMinLat = loadLat# 更新最大纬度if loadLat > resultMaxLat:resultMaxLat = loadLat# 结果输出
print("厦门最小经度：", resultMinLon)
print("厦门最大经度：", resultMaxLon)
print("厦门最小纬度：", resultMinLat)
print("厦门最大纬度：", resultMaxLat)# 可视化
fig, ax = plt.subplots(figsize=(10, 10))# 绘制厦门市的边界
loadGeoData.plot(ax=ax, edgecolor='black', facecolor='lightblue')# 绘制最小外接矩形
min_rect = Rectangle((resultMinLon, resultMinLat), resultMaxLon - resultMinLon, resultMaxLat - resultMinLat,linewidth=2, edgecolor='red', facecolor='none', linestyle='--')
ax.add_patch(min_rect)# 标注最小和最大经度与纬度点
ax.scatter(resultMinLon, resultMinLat, color='red', label='最小边界')
ax.scatter(resultMaxLon, resultMaxLat, color='green', label='最大边界')
ax.annotate(f'({resultMinLon:.6f}, {resultMinLat:.6f})', (resultMinLon, resultMinLat), textcoords="offset points",xytext=(0, 10), ha='center')
ax.annotate(f'({resultMaxLon:.6f}, {resultMaxLat:.6f})', (resultMaxLon, resultMaxLat), textcoords="offset points",xytext=(0, 10), ha='center')# 添加图例
ax.legend()# 设置坐标轴范围
ax.set_xlim(resultMinLon - 0.01, resultMaxLon + 0.01)
ax.set_ylim(resultMinLat - 0.01, resultMaxLat + 0.01)# 显示图形
plt.title('厦门市边界及最小外接矩形')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.grid(True)
plt.show()

下面这部分代码增加了一个功能：就是把这个矩形导出来，这里是把这个矩形导出来形成shp文件，修改文件路径的基础上再改一下保存位置路径即可；

完整代码#运行环境Python 3.11

python">import geopandas as gp
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from shapely.geometry import box# 配置 matplotlib 支持中文显示
plt.rcParams['font.family'] = ['SimHei']  # 使用黑体字体
plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号# 结果变量
resultMinLon = 99999.0  # 最小经度
resultMaxLon = -99999.0  # 最大经度
resultMinLat = 99999.0  # 最小纬度
resultMaxLat = -99999.0  # 最大纬度# 读取厦门市边界数据
# 注意：请确保文件路径正确
loadGeoData = gp.read_file("D:\\data\\厦门市.json")# 遍历每一行数据
for i in range(len(loadGeoData)):# 读取几何数据loadGeometry = loadGeoData.loc[i, "geometry"]# 检查几何数据类型，可能是单个多边形（Polygon）或多个多边形（MultiPolygon）if loadGeometry.geom_type == 'Polygon':# 单个多边形# 遍历多边形的所有顶点for z in range(len(loadGeometry.exterior.coords)):# 获取当前顶点的经度和纬度loadLon, loadLat = loadGeometry.exterior.coords[z]# 更新最小经度if loadLon < resultMinLon:resultMinLon = loadLon# 更新最大经度if loadLon > resultMaxLon:resultMaxLon = loadLon# 更新最小纬度if loadLat < resultMinLat:resultMinLat = loadLat# 更新最大纬度if loadLat > resultMaxLat:resultMaxLat = loadLatelif loadGeometry.geom_type == 'MultiPolygon':# 多个多边形# 遍历每个多边形for polygon in loadGeometry.geoms:# 遍历多边形的所有顶点for z in range(len(polygon.exterior.coords)):# 获取当前顶点的经度和纬度loadLon, loadLat = polygon.exterior.coords[z]# 更新最小经度if loadLon < resultMinLon:resultMinLon = loadLon# 更新最大经度if loadLon > resultMaxLon:resultMaxLon = loadLon# 更新最小纬度if loadLat < resultMinLat:resultMinLat = loadLat# 更新最大纬度if loadLat > resultMaxLat:resultMaxLat = loadLat# 结果输出
print("厦门最小经度：", resultMinLon)
print("厦门最大经度：", resultMaxLon)
print("厦门最小纬度：", resultMinLat)
print("厦门最大纬度：", resultMaxLat)# 创建最小外接矩形的几何对象
bounding_box = box(resultMinLon, resultMinLat, resultMaxLon, resultMaxLat)# 创建 GeoDataFrame
bbox_gdf = gp.GeoDataFrame(geometry=[bounding_box], crs=loadGeoData.crs)# 保存为 Shapefile
output_path = "D:\\data\\厦门最小外接矩形.shp"
bbox_gdf.to_file(output_path)# 可视化
fig, ax = plt.subplots(figsize=(10, 10))# 绘制厦门市的边界
loadGeoData.plot(ax=ax, edgecolor='black', facecolor='lightblue')# 绘制最小外接矩形
min_rect = Rectangle((resultMinLon, resultMinLat), resultMaxLon - resultMinLon, resultMaxLat - resultMinLat,linewidth=2, edgecolor='red', facecolor='none', linestyle='--')
ax.add_patch(min_rect)# 标注最小和最大经度与纬度点
ax.scatter(resultMinLon, resultMinLat, color='red', label='最小边界')
ax.scatter(resultMaxLon, resultMaxLat, color='green', label='最大边界')
ax.annotate(f'({resultMinLon:.6f}, {resultMinLat:.6f})', (resultMinLon, resultMinLat), textcoords="offset points",xytext=(0, 10), ha='center')
ax.annotate(f'({resultMaxLon:.6f}, {resultMaxLat:.6f})', (resultMaxLon, resultMaxLat), textcoords="offset points",xytext=(0, 10), ha='center')# 添加图例
ax.legend()# 设置坐标轴范围
ax.set_xlim(resultMinLon - 0.01, resultMaxLon + 0.01)
ax.set_ylim(resultMinLat - 0.01, resultMaxLat + 0.01)# 显示图形
plt.title('厦门市边界及最小外接矩形')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.grid(True)
plt.show()

忽然发现arcgis/arcgispro也可以生成这个最小外接矩形，直接检索【最小边界几何】，几何类型选择【按面积矩形】或者【按宽度矩形】都可，看需求，组选型选【ALL】这样会基于全部结果形成最外层的矩形；