在 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. To get some inspiration, you can also verify the matlab function areaintn 的代码。我不知道它是否直接存在于 python.

中的类似函数

Answer 3

我正在利用给出的答案寻找两个旋转矩形的交点（请找到原始答案）。如果您觉得这个答案有用，请不要给它点赞，但请去给原始海报点赞。我不认为这是我的功劳。此外，此答案必须适应您的具体情况。

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

基于其他人发布的所有 answers/comments，最终当我在末尾添加这些行时它起作用了。没有回复者的帮助，我做不到：

#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,'%')

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

Calculate overlap between polygon and shapefile in Python 3.6

python

cartopy