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

Published: October 30, 2023

## 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`
```

Image

of