Introduction
In this article, we will learn how to find the intersection between two geographic areas using two popular python libraries - shapely and cartopy.
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:
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 |