Calculate overlap between polygon and shapefile in Python 3.6

Question

I would like to calculate the percentage of overlap between a shapefile and a polygon. I'm using Cartopy and Matplotlib and created the map shown here:

A part of Europe (using a shapefile downloaded here) and an arbitrary rectangle are shown. Let's say I would like to calculate the percentage of Belgium that is covered by the rectangle. How would I do this? Below, the code so far is shown.

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

score 8 · Accepted Answer · answered May 18 '18 at 08:11

8

Based on all the answers/comments posted by others, eventually it worked when I added these lines at the end. I could not have done it without the help of the respondees:

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

answered May 18 '18 at 08:11

A T

281
3
11

1

Nice you got your answer. Remember that "variable.area" will tell you the area in the euclidian space, and what you're doing it's a polygon on the surface of the sphere, which might considerably change the result. – silgon May 21 '18 at 14:26
Is there a way to get the area for a curved (non-Euclidian) surface? – A T May 22 '18 at 07:53

silgon · Answer 2 · 2018-05-16T14:13:16.103

5

If you have the shape of both polygons. You can do something like the following:

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

Of course, I'm oversimplifying the problem and the result most probably is with euclidian area which might not be what you need. Thus said, for the area of a polygon on the surface of a sphere you can verify code from the following reference. To get some inspiration, you can also verify the matlab function areaintn. I don't know if it exists directly an analogous function in python.

edited May 16 '18 at 14:13

answered May 16 '18 at 14:06

silgon

6,890
7
46
67

The problem is that the shapefile is not a polygon (I think). When using `PLT.intersection(poly)`, it yields the following error: `AttributeError: 'GeoAxesSubplot' object has no attribute 'intersection'` – A T May 16 '18 at 14:15
1

@AT the trivial solution here is that you have to find a way to "close" your polygon in ArcGIS or QGIS or whatever you are using, or find a Pythonic solution to do that. – FaCoffee May 16 '18 at 14:16
1

@AT you can transform it, look at [this reference](http://scitools.org.uk/cartopy/docs/v0.15/tutorials/using_the_shapereader.html) and [this reference](http://scitools.org.uk/cartopy/docs/v0.14/crs/index.html). Look for shapely in both references, it's kind of standard. – silgon May 16 '18 at 14:18

score 2 · Answer 3 · answered May 16 '18 at 13:28

2

Well, you'd need some sort of way of getting the area of the entire map, as well as the area of the country itself. The area of the polygon is probably the most simple part.

I'd suggest start with something more basic, maybe just a grid and two simple shapes, this could help you envision how it could be done at a more complex level.

answered May 16 '18 at 13:28

Chandler

21
2

Indeed, the area of the polygon is simply found by `print (poly.area)`. I do agree that this is an advanced problem, but it would already greatly help me to find the area of e.g. Belgium. – A T May 16 '18 at 13:38

score 1 · Answer 4 · answered May 16 '18 at 14:14

I am taking advantage of the answer given for finding the intersection of two rotated rectangles (please find the original answer here). Please don't upvote this answer if you find it useful, but go and upvote the original poster. I do not take credit for this. Also, this answer must be adapted to your specific case.

TL:DR The answer involves using 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()

Calculate overlap between polygon and shapefile in Python 3.6

4 Answers4