如何在Python 3.6中计算多边形与Shapefile的重叠占比？

嘿，要算出比利时被你的自定义矩形覆盖的面积占比，咱们得调整一下现有代码，重点解决两个问题：准确提取比利时的边界，以及用合适的投影计算面积——毕竟经纬度坐标系下算面积会有误差。下面是具体的实现方案：

核心思路 #

1. 从Shapefile中精准筛选出比利时的几何图形（不用加载所有国家）
2. 把比利时边界和自定义矩形都转换到等面积投影（比如欧洲常用的EPSG:3035），确保面积计算准确
3. 计算比利时总面积、重叠区域面积，最后得出占比

完整代码 #

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shapereader
from shapely.geometry import Polygon
from descartes import PolygonPatch
from shapely.ops import transform
from pyproj import Transformer

# 1. 初始化绘图设置
fig1 = plt.figure(figsize=(10,10))
PLT = plt.axes(projection=ccrs.PlateCarree())
PLT.set_extent([-10,10,45,55])
PLT.gridlines()

# 2. 加载Shapefile并提取比利时的几何边界
fname = r'C:\Users\Me\ne_50m_admin_0_countries.shp'
reader = shapereader.Reader(fname)
belgium_geom = None

# 遍历所有国家记录，匹配比利时的名称
for record in reader.records():
    if record.attributes['NAME'] == 'Belgium':
        belgium_geom = record.geometry
        break

# 绘制所有国家基础图层，并高亮比利时
adm1_shapes = list(reader.geometries())
PLT.add_geometries(adm1_shapes, ccrs.PlateCarree(), edgecolor='black', facecolor='gray', alpha=0.5)
if belgium_geom:
    PLT.add_geometries([belgium_geom], ccrs.PlateCarree(), edgecolor='red', facecolor='orange', alpha=0.7)

# 3. 创建自定义矩形多边形
x3 = 4
x4 = 5
y3 = 50
y4 = 52
poly = Polygon([(x3,y3),(x3,y4),(x4,y4),(x4,y3)])
PLT.add_patch(PolygonPatch(poly, fc='#cc00cc', ec='#555555', alpha=0.5, zorder=5))

# 4. 投影转换：用等面积投影计算面积（避免经纬度的面积误差）
# 转换器：从WGS84经纬度（EPSG:4326）转欧洲等面积投影（EPSG:3035）
transformer = Transformer.from_crs("EPSG:4326", "EPSG:3035", always_xy=True)

# 转换两个几何对象到等面积投影
belgium_equal_area = transform(transformer.transform, belgium_geom)
poly_equal_area = transform(transformer.transform, poly)

# 5. 计算面积占比
if belgium_equal_area and poly_equal_area:
    belgium_total_area = belgium_equal_area.area
    overlap_area = belgium_equal_area.intersection(poly_equal_area).area
    overlap_ratio = (overlap_area / belgium_total_area) * 100

    # 打印结果（转换为平方公里更易读）
    print(f"比利时被矩形覆盖的面积占比：{overlap_ratio:.2f}%")
    print(f"比利时总面积：{belgium_total_area / 1e6:.2f} 平方公里")
    print(f"重叠区域面积：{overlap_area / 1e6:.2f} 平方公里")

plt.show()

关键细节说明 #

• 投影转换的必要性：咱们用的ccrs.PlateCarree()是经纬度坐标系，直接算面积会因为高纬度的经度拉伸导致结果不准。EPSG:3035是专门为欧洲设计的等面积投影，能保证面积计算的准确性。
• 提取比利时边界：通过遍历Shapefile的属性记录，匹配NAME字段为"Belgium"的条目，这样就能精准获取比利时的几何形状。
• 重叠区域计算：用Shapely的intersection()方法获取两个几何的重叠部分，再通过area属性计算面积，最后算出占比。

内容的提问来源于stack exchange，提问作者A T