在 Python 3.6 中计算多边形和 shapefile 之间的重叠

Question

我想计算 shapefile 和多边形之间的重叠百分比。我正在使用 Cartopy 和 Matplotlib 并创建了此处显示的地图：

显示了欧洲的一部分（使用下载的 shapefile here）和任意矩形。假设我想计算矩形覆盖的比利时的百分比。我该怎么做？下面是目前为止的代码。

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

#create figure
fig1 = plt.figure(figsize=(10,10))  
PLT = plt.axes(projection=ccrs.PlateCarree())
PLT.set_extent([-10,10,45,55])
PLT.gridlines()

#import and display shapefile
fname = r'C:\Users\Me\ne_50m_admin_0_countries.shp'
adm1_shapes = list(shapereader.Reader(fname).geometries())
PLT.add_geometries(adm1_shapes, ccrs.PlateCarree(),
              edgecolor='black', facecolor='gray', alpha=0.5)

#create arbitrary polygon
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))

Answer 1

嗯，你需要某种方法来获取整个地图的面积，以及国家本身的面积。多边形的面积可能是最简单的部分。

我建议从更基本的东西开始，也许只是一个网格和两个简单​​的形状，这可以帮助您设想如何在更复杂的层面上完成。

Answer 2

如果你有两个多边形的形状。您可以执行以下操作：

from shapely.geometry import Polygon
p1=Polygon([(0,0),(1,1),(1,0)])
p2=Polygon([(0,1),(1,0),(1,1)])
p3=p1.intersection(p2)
print(p3) # result: POLYGON ((0.5 0.5, 1 1, 1 0, 0.5 0.5))
print(p3.area) # result: 0.25

当然，我把问题简单化了，结果很可能是欧几里得面积，这可能不是你需要的。因此，对于球体表面上的多边形区域，您可以验证来自following reference 的代码。为了得到一些启发，还可以验证matlab函数areaintn。不知道是不是直接存在python中的类似函数。

Answer 3

我正在利用给出的答案来寻找两个旋转矩形的交点（请找到原始答案here）。如果您觉得它有用，请不要赞成这个答案，而是去赞成原始海报。我不相信这一点。此外，此答案必须适合您的具体情况。

TL:DR 答案涉及使用shapely。

    import shapely.geometry
    import shapely.affinity

    class RotatedRect:
        def __init__(self, cx, cy, w, h, angle):
            self.cx = cx
            self.cy = cy
            self.w = w
            self.h = h
            self.angle = angle

        def get_contour(self):
            w = self.w
            h = self.h
            c = shapely.geometry.box(-w/2.0, -h/2.0, w/2.0, h/2.0)
            rc = shapely.affinity.rotate(c, self.angle)
            return shapely.affinity.translate(rc, self.cx, self.cy)

        def intersection(self, other):
            return self.get_contour().intersection(other.get_contour())


    r1 = RotatedRect(10, 15, 15, 10, 30)
    r2 = RotatedRect(15, 15, 20, 10, 0)

    from matplotlib import pyplot
    from descartes import PolygonPatch

    fig = pyplot.figure(1, figsize=(10, 4))
    ax = fig.add_subplot(121)
    ax.set_xlim(0, 30)
    ax.set_ylim(0, 30)

    ax.add_patch(PolygonPatch(r1.get_contour(), fc='#990000', alpha=0.7))
    ax.add_patch(PolygonPatch(r2.get_contour(), fc='#000099', alpha=0.7))
    ax.add_patch(PolygonPatch(r1.intersection(r2), fc='#009900', alpha=1))

    pyplot.show()

Answer 4

根据其他人发布的所有答案/cmets，当我在最后添加这些行时，它最终起作用了。如果没有受访者的帮助，我是做不到的：

#find the Belgium polygon. 
for country in shapereader.Reader(fname).records(): 
    if country.attributes['SOV_A3'] == 'BEL': 
        Belgium = country.geometry 
        break 
PLT.add_patch(PolygonPatch(poly.intersection(Belgium), fc='#009900', alpha=1)) 

#calculate coverage 
x = poly.intersection(Belgium) 
print ('Coverage: ', x.area/Belgium.area*100,'%')