How to check if two geographical areas intersect in python using shapely ?

Introduction

In this document, we will be discussing how to check if two geographical areas intersect in Python using the shapely library. Before diving into the technical details, let's first understand what is meant by geographical areas and why it is important to check for their intersections.

Understanding Geographical Areas

Geographical areas refer to regions on the Earth's surface that are defined by their boundaries, such as countries, cities, states, or even smaller subdivisions like neighborhoods. They can also be represented by polygons on a map.

Checking Intersections with Shapely

Shapely is a popular Python library used for geometric operations on Cartesian plane figures. It provides functions to create and manipulate geometric objects like points, lines, and polygons.

To check if two geographical areas intersect using shapely, we first need to create two polygon objects representing the respective areas. Then, we can use the intersects() function to determine if they intersect with each other. If the function returns a True value, it means that there is an intersection between the two polygons.

Import necessary libraries

from shapely.geometry import Point
from shapely.ops import transform

import pyproj

Create two circles

pt1 = Point(2.3522, 48.85).buffer(10) # longitude, latitde, radius
pt2 = Point(37.6173, 55.7558).buffer(20)

Convert Cartesian coordinates to a geodetic projection

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

pt1 = transform(project.transform, pt1)
pt2 = transform(project.transform, pt2)

Checking Intersections

pt1.intersects(pt2)

The output of this code would be

False

Another example

from shapely.geometry import Point
from shapely.ops import transform

import pyproj

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

pt1 = Point(2.3522, 48.85).buffer(15) # longitude, latitde
pt2 = Point(37.6173, 55.7558).buffer(25)

pt1 = transform(project.transform, pt1)
pt2 = transform(project.transform, pt2)

pt1.intersects(pt2)

The output of this code would be

True

Visualizing areas of intersection using cartopy

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='110m')

poly = mpatches.Circle((2.3522, 48.85), 10.0, color='coral', transform=ccrs.PlateCarree())
ax.add_patch(poly)

poly = mpatches.Circle((37.6173, 55.7558), 20.0, color='lightblue', transform=ccrs.PlateCarree())
ax.add_patch(poly)

plt.title('How to check if two geographical areas intersect in python using shapely ?')

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

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='110m')

poly = mpatches.Circle((2.3522, 48.85), 15.0, color='coral', transform=ccrs.PlateCarree())
ax.add_patch(poly)

poly = mpatches.Circle((37.6173, 55.7558), 25.0, color='lightblue', transform=ccrs.PlateCarree(), alpha=0.5)
ax.add_patch(poly)

plt.title('How to check if two geographical areas intersect in python using shapely ?')

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

Another example: Intersection between two polygons

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

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', fill=True, lw=1, fc=None, transform=ccrs.Geodetic(), alpha=0.5)
ax.add_patch(poly)

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', fill=True, lw=1, fc=None, transform=ccrs.Geodetic(), alpha=0.5)
ax.add_patch(poly)

plt.title('How to check if two geographical areas intersect in python using shapely ?')

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

plt.show()

from shapely.geometry import Polygon

import pyproj

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

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

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

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

polygon2 = Polygon(lons_lats_vect) # create polygon

polygon1.intersects(polygon2)

The output of this code would be

True

References