当前位置： 首页 > article >正文

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

article 2024/11/13 16:08:31

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

结果如下；

文章仅用于分享个人学习成果与个人存档之用，分享知识，如有侵权，请联系作者进行删除。所有信息均基于作者的个人理解和经验，不代表任何官方立场或权威解读。


http://www.kler.cn/a/306662.html

相关文章：

  • 怎么监控员工电脑？分享5个监控员工电脑的绝佳方法（立竿见影！建议收藏！）
  • 【Threejs】相机控制器动画
  • ELK-Logstash配置
  • vscode Markdown
  • 前端入门一之ES6--面向对象、够着函数和原型、继承、ES5新增方法、函数进阶、严格模式、高阶函数、闭包
  • 希尔排序（C语言）
  • ubuntu 安装 chrome 及 版本匹配的 chromedriver
  • vue3+vite项目中使用阿里图标库(svg)图标
  • NX CAM二次开发-创建程序组
  • Linux套接字
  • Python Web 开发中的性能优化策略（一）
  • Java多线程面试精讲：源于技术书籍的深度解读
  • uniapp+vue3 使用canvas，并保存图片（图片是空白的问题）
  • PMP–一、二、三模–分类–14.敏捷–技巧–项目生命周期
  • LINUX网络编程:http
  • HSmartWindowControl 滚轮缩放 交互式绘制ROI 可修改 存储
  • 初写MySQL四张表:(2/4)
  • 切换淘宝最新npm镜像源
  • Android应用性能优化
  • 抚琴成一快-音程和度数
  • 证券api接口，一个开源Python量化交易平台项目需要考虑哪些方面
  • [JVM]JVM内存划分, 类加载过程, 双亲委派模型,垃圾回收机制
  • 学习笔记JVM篇（一）
  • C语言中的信号量应用
  • 【ArcGIS Pro实操第七期】栅格数据合并、裁剪及统计：以全球不透水面积为例
  • Linux03