How to find the intersection between two geographic areas in python using shapely and cartopy ?

Introduction

In this article, we will learn how to find the intersection between two geographic areas using two popular python libraries - shapely and cartopy.

Table of contents

Introduction
Creating two polygons using shapely
Transforming the coordinates into a projected coordinate system
Extracting the intersection
Inverse transformation
Visualizations
- Plotting the intersection
- Plotting all together
References

Shapely is a python library used for manipulating and analyzing geometric objects. It provides various functions to perform operations such as intersection, union, difference, and more on geometric shapes. On the other hand, cartopy is a mapping library built upon matplotlib that provides a simple and powerful way to handle geospatial data. It allows us to create maps, plot spatial data, and perform various transformations. Together, these two libraries provide a robust solution for finding the intersection between geographic areas.

Let's explore a real-life case study where we are presented with two images from a NOAA satellite. Our objective is to identify and determine the overlapping area between these images.

To gain a deeper understanding of our problem, let's begin by visualizing our dataset:

from matplotlib.pyplot import figure

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib.patches as mpatches

fig = figure(num=None, figsize=(12, 10), dpi=100, edgecolor='k')

map_proj = ccrs.Robinson()

ax = plt.axes(projection=map_proj)

ax.set_global()
ax.gridlines()

ax.coastlines(linewidth=0.5, color='k', resolution='50m')

ax.set_extent([-16.0, 52.0, 43.0, 55.5])

#-----------------------------------------------------------------#
# Satellite Obs 1

lat_corners = np.array([47.956543, 54.523647, 49.484398, 43.507904])
lon_corners = np.array([52.530487, 9.432528, 9.602138, 48.74224])

poly_corners = np.zeros((len(lat_corners), 2), np.float64)
poly_corners[:,0] = lon_corners
poly_corners[:,1] = lat_corners

poly = mpatches.Polygon(poly_corners, 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='coral', 
                        transform=ccrs.Geodetic(), 
                        alpha=0.5)

ax.add_patch(poly)

#-----------------------------------------------------------------#
# Satellite Obs 2

lat_corners = np.array([49.427364, 56.230656, 51.196304, 45.03514])
lon_corners = np.array([28.614933, -16.069532, -15.793909, 24.584558])

poly_corners = np.zeros((len(lat_corners), 2), np.float64)
poly_corners[:,0] = lon_corners
poly_corners[:,1] = lat_corners

poly = mpatches.Polygon(poly_corners, 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='coral',
                        transform=ccrs.Geodetic(), 
                        alpha=0.5)

ax.add_patch(poly)

#-----------------------------------------------------------------#

plt.title('How to find the intersection between two geographic areas in python using shapely and cartopy ?')

plt.savefig("viirs_granules_intersection_01.png", bbox_inches='tight', facecolor='white')

plt.show()

The above code will generate the following image:

How to find the intersection between two geographic areas in python using shapely and cartopy ?

Creating two polygons using shapely

First, we can represent the geographic boundaries of each image using shapely.

Image 1:

from shapely.geometry import Polygon

lat_corners = np.array([47.956543, 54.523647, 49.484398, 43.507904])
lon_corners = np.array([52.530487, 9.432528, 9.602138, 48.74224])

lons_lats_vect = np.column_stack((lon_corners, lat_corners)) # Reshape coordinates

polygon1 = Polygon(lons_lats_vect) # create polygon

Image 2:

lat_corners = np.array([49.427364, 56.230656, 51.196304, 45.03514])
lon_corners = np.array([28.614933, -16.069532, -15.793909, 24.584558])

#lat_corners = np.array([25.303867, 29.96691, 24.902569, 20.422632])
#lon_corners = np.array([38.612545, 8.326748, 7.698123, 36.761963])

lons_lats_vect = np.column_stack((lon_corners, lat_corners)) # Reshape coordinates

polygon2 = Polygon(lons_lats_vect) # create polygon

Transforming the coordinates into a projected coordinate system

In order to determine the intersection between the two images, it is crucial to convert the coordinates of the images, namely longitude and latitude, into a projected coordinate system. This particular step holds great significance as the calculation of the intersection can only be performed using planar coordinates:

from shapely.ops import transform

import pyproj

project = pyproj.Transformer.from_proj(
    pyproj.Proj("EPSG:4326"), 
    pyproj.Proj('EPSG:25832'),
    always_xy=True
)

polygon1 = transform(project.transform, polygon1)
polygon2 = transform(project.transform, polygon2)

Extracting the intersection

To find the intersection between two geometric objects, we will use the intersection() function from shapely. This function returns a new geometric object that represents the overlapping region between the two inputs. In our case, we will pass in our county and park shapes as follows:

polygon3 = polygon1.intersection(polygon2)

To obtain the coordinates of the new polygon, one possible solution is to utilize the exterior.coords.xy method

X,Y = polygon3.exterior.coords.xy

X,Y

gives

(array('d', [1880230.956486944, 543613.0916985297, 527994.9156186895, 591345.3283386811, 1917597.861425543, 1880230.956486944]),
 array('d', [5553405.894095841, 5481480.503905152, 6041871.479147676, 6047056.38475262, 5662676.336762555, 5553405.894095841]))

Inverse transformation

To obtain the latitudes and longitudes of the intersection, we can simply apply the inverse transformation

project = pyproj.Transformer.from_proj(
    pyproj.Proj('EPSG:25832'),
    pyproj.Proj("EPSG:4326"), 
    always_xy=True
)

polygon4 = transform(project.transform, polygon3)

Then

X,Y = polygon4.exterior.coords.xy

X,Y

gives

(array('d', [27.761722057952035, 9.602138, 9.432528, 10.412699947238812, 28.61493299999999, 27.761722057952035]),
 array('d', [48.58246264585949, 49.48439799999999, 54.52364699999999, 54.56277121439247, 49.42736399999999, 48.58246264585949]))

Which can be used to visualize the result.

Visualizations

Plotting the intersection

from matplotlib.pyplot import figure

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib.patches as mpatches

fig = figure(num=None, figsize=(12, 10), dpi=100, edgecolor='k')

map_proj = ccrs.Robinson()

ax = plt.axes(projection=map_proj)

ax.set_global()
ax.gridlines()

ax.coastlines(linewidth=0.5, color='k', resolution='50m')

ax.set_extent([-16.0, 52.0, 43.0, 55.5])

lat_corners = np.array([48.58246264585949, 49.48439799999999, 54.52364699999999, 54.56277121439247, 49.42736399999999, 48.58246264585949])
lon_corners = np.array([27.761722057952035, 9.602138, 9.432528, 10.412699947238812, 28.61493299999999, 27.761722057952035])

poly_corners = np.zeros((len(lat_corners), 2), np.float64)
poly_corners[:,0] = lon_corners
poly_corners[:,1] = lat_corners

poly = mpatches.Polygon(poly_corners, 
                        closed=False, 
                        ec='k', 
                        lw=0.4, 
                        color='darkgray',
                        transform=ccrs.Geodetic())

ax.add_patch(poly)

#-----------------------------------------------------------------#

plt.title('How to find the intersection between two geographic areas in python using shapely and cartopy ?')

plt.savefig("viirs_granules_intersection_02.png", bbox_inches='tight', facecolor='white')

plt.show()

Plotting all together

from matplotlib.pyplot import figure

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib.patches as mpatches

fig = figure(num=None, figsize=(12, 10), dpi=100, edgecolor='k')

map_proj = ccrs.Robinson()

ax = plt.axes(projection=map_proj)

ax.set_global()
ax.gridlines()

ax.coastlines(linewidth=0.5, color='k', resolution='50m')

ax.set_extent([-16.0, 52.0, 43.0, 55.5])

#-----------------------------------------------------------------#
# Satellite Obs 1

lat_corners = np.array([47.956543, 54.523647, 49.484398, 43.507904])
lon_corners = np.array([52.530487, 9.432528, 9.602138, 48.74224])

poly_corners = np.zeros((len(lat_corners), 2), np.float64)
poly_corners[:,0] = lon_corners
poly_corners[:,1] = lat_corners

poly = mpatches.Polygon(poly_corners, 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='coral', 
                        transform=ccrs.Geodetic(), 
                        alpha=0.5)

ax.add_patch(poly)

#-----------------------------------------------------------------#
# Satellite Obs 2

lat_corners = np.array([49.427364, 56.230656, 51.196304, 45.03514])
lon_corners = np.array([28.614933, -16.069532, -15.793909, 24.584558])

poly_corners = np.zeros((len(lat_corners), 2), np.float64)
poly_corners[:,0] = lon_corners
poly_corners[:,1] = lat_corners

poly = mpatches.Polygon(poly_corners, 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='coral',
                        transform=ccrs.Geodetic(), 
                        alpha=0.5)

ax.add_patch(poly)

#-----------------------------------------------------------------#
# Intersection

lat_corners = np.array([48.58246264585949, 49.48439799999999, 54.52364699999999, 54.56277121439247, 49.42736399999999, 48.58246264585949])
lon_corners = np.array([27.761722057952035, 9.602138, 9.432528, 10.412699947238812, 28.61493299999999, 27.761722057952035])

poly_corners = np.zeros((len(lat_corners), 2), np.float64)
poly_corners[:,0] = lon_corners
poly_corners[:,1] = lat_corners

poly = mpatches.Polygon(poly_corners, 
                        closed=False, 
                        ec='k', 
                        lw=0.4, 
                        color='darkgray',
                        transform=ccrs.Geodetic())

ax.add_patch(poly)

#-----------------------------------------------------------------#

plt.title('How to find the intersection between two geographic areas in python using shapely and cartopy ?')

plt.savefig("viirs_granules_intersection_03.png", bbox_inches='tight', facecolor='white')

plt.show()

References

Links	Site
How to check if two geographical areas intersect in python using shapely ?	moonbooks.org
How to plot a circle on a cartopy map with python ?	moonbooks.org
How to extract the intersection between two polygons in Python using Shapely ?	moonbooks.org

Developer in my spare time, I’m interested in Earth observation and satellite-based fire detection. Passionate about Python, machine learning, and open science.