Вычислить перекрытие между полигоном и шейп-файлом в 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 areaint n. Я не знаю, существует ли в 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

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

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