有的时候我们用行政区的面层进行裁剪的时候，特别是在基于行政区裁剪的渔网，因为行政边界本身就是不规则的原因，在边界的渔网匹配数据会出现匹配不上的问题，其实如果周边有数据完全可以采用最近匹配的原则，但是就是因为行政区裁太可丁可卯了，会丢失一部分数据，这个时候可以通过形成一个最小外接矩形的形式，可以让数据裁剪保持一个微妙的范围，不至于太大，又不至于没有二次处理的空间，第二个方面就是比如我们需要行政区里面的POI，以高德为例，因为行政区本身不是规则的，我们通过建立外接矩阵，这样获取其中的数据就不会出现丢失的情况了；
行政区边界获取部分，这里采用了阿里云旗下的数据可视化平台，其数据源本身也是来着高德地图，数据质量还是比较可靠的；
全国行政区划边界GeoJSON数据下载平台阿里云：DataV.GeoAtlas地理小工具系列 (aliyun.com)
这里把下载到本地的geojson或者json的文件路径改成自己的即可；
完整代码#运行环境Python 3.11
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 = loadLat
elif 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
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 = loadLat
elif 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
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 = loadLat
elif 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】这样会基于全部结果形成最外层的矩形；
结果如下；
文章仅用于分享个人学习成果与个人存档之用，分享知识，如有侵权，请联系作者进行删除。所有信息均基于作者的个人理解和经验，不代表任何官方立场或权威解读。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
