Вычислить перекрытие между полигоном и шейп-файлом в Python 3.6

Question

Я хотел бы рассчитать процент перекрытия между шейп-файлом и полигоном. Я использую Cartopy и Matplotlib и создал карту, показанную здесь:

Часть Европы (с использованием шейп-файла, загруженного здесь) и произвольный прямоугольник показаны. Допустим, я хотел бы рассчитать процент Бельгии, который покрыт прямоугольником. Как бы я это сделал? Ниже пока показан код.

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

Основываясь на всех ответах / комментариях, опубликованных другими, в конце концов это сработало, когда я добавил эти строки в конце. Я не мог бы сделать это без помощи респондентов:

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

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

Конечно, я упрощаю проблему, и, скорее всего, результат будет с евклидовой областью, которая может оказаться не той, что вам нужна. Таким образом, для области многоугольника на поверхности сферы вы можете проверить код из следующей ссылки. Чтобы получить вдохновение, вы также можете проверить, что функция Matlab незнакома. Я не знаю, существует ли напрямую аналогичная функция в Python.

Answer 3

Ну, вам нужен какой-то способ получить площадь всей карты, а также площадь самой страны. Площадь многоугольника, вероятно, самая простая часть.

Я бы предложил начать с чего-то более простого, может быть, просто с сетки и двух простых фигур, это может помочь вам представить, как это можно сделать на более сложном уровне.

Answer 4

Я использую ответ, данный для нахождения пересечения двух повернутых прямоугольников (пожалуйста, найдите оригинальный ответ здесь). Пожалуйста, не голосуйте за этот ответ, если вы считаете его полезным, но идите и проголосуйте за оригинальный постер. Я не беру кредит на это. Кроме того, этот ответ должен быть адаптирован к вашему конкретному случаю.

TL: DR. Ответ предполагает использование стройного.

    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()