How to determine if a point, specified by its latitude and longitude, falls within a given area using the shapely library in Python ?

Using the Shapely library in Python, it is possible to check if a point with a latitude and longitude is within an area or not. First, you need to create a shapely Geometry object for your area of interest (polygon). Then you can create a shapely Point object with your latitude and longitude. Finally, you can check if this point is within the area defined by the polygon using the contains or within methods, examples:

Installing Shapely

Shapely can be installed, for example, along with geopandas as it is a dependency:

pip install geopandas

Check if a point is within a given polygon

Please import the required Python libraries:

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

import numpy as np

Define a shapely point

Define a target point, for example

target_lat = 43.7580
target_lon = -65.3208

point = Point(target_lon,target_lat) # create point

Define a shapely polygon

Define an area using a shapely polygon

lons_vect = np.array([-83.57072,-47.520203,-47.604248,-86.768524])
lats_vect = np.array([38.963406,44.391663,49.442635,43.54688,])

lons_lats_vect = np.column_stack((lons_vect, lats_vect)) # Reshape coordinates

polygon = Polygon(lons_lats_vect) # create polygon

Using shapely contains() method

Now you can use this method to check if a point with a latitude and longitude is within an area of interest. It is important to note that the coordinates should be provided in the same coordinate system as used for your polygon or else it will produce erroneous results.:

polygon.contains(point)

returns here

True

Using shapely within() method

To check if a point is in the polygon

point.within(polygon)

returns here

True

Data visualization

import matplotlib.pyplot as plt
import matplotlib as mpl
import cartopy.crs as ccrs
import matplotlib.patches as mpatches

plt.figure(figsize=(16,9))

proj = ccrs.PlateCarree()

ease_extent = [-180., 180., 90., -90.]

ax = plt.axes(projection=proj)

ax.set_extent(ease_extent, crs=proj)

ax.gridlines(color='gray', linestyle='--')

ax.coastlines()

lat_corners = lats_vect
lon_corners = lons_vect

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, fill=False, lw=1, fc=None, transform=ccrs.PlateCarree())
ax.add_patch(poly)

plt.scatter( target_lon, target_lat,
         color='red', linewidth=2, marker='o', s=20,
         transform=ccrs.PlateCarree(),
         )

plt.scatter( lons_vect, lats_vect,
         color='blue', linewidth=2, marker='o', s=20,
         transform=ccrs.PlateCarree(),
         )

plt.tight_layout()

plt.show()

plt.close()

Transforming Shapely Polygon and point in another coordinate system

Now let's consider the following example

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

import numpy as np

target_lat =  47.18241500854492
target_lon = -115.90896606445312

lons_vect = np.array([-104.50501251, -140.47285461, -143.56848145, -104.57888031])
lats_vect = np.array([43.46247482, 38.01778412, 42.60084152, 48.49499893])

lons_lats_vect = np.column_stack((lons_vect, lats_vect)) # Reshape coordinates

polygon = Polygon(lons_lats_vect) # create polygon

point = Point(target_lon,target_lat) # create point

print(polygon.contains(point)) # check if polygon contains point
print(point.within(polygon)) # check if a point is in the polygon

it will returns

False
False

Let's plot the data

import matplotlib.pyplot as plt
import matplotlib as mpl
import cartopy.crs as ccrs
import matplotlib.patches as mpatches

plt.figure(figsize=(16,9))

proj = ccrs.PlateCarree()

ease_extent = [-180., 180., 90., -90.]

ax = plt.axes(projection=proj)

ax.set_extent(ease_extent, crs=proj)

ax.gridlines(color='gray', linestyle='--')

ax.coastlines()

lat_corners = lats_vect
lon_corners = lons_vect

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, fill=False, lw=1, fc=None, transform=ccrs.PlateCarree())
ax.add_patch(poly)

plt.scatter( target_lon, target_lat,
         color='red', linewidth=2, marker='o', s=20,
         transform=ccrs.PlateCarree(),
         )

plt.scatter( lons_vect, lats_vect,
         color='blue', linewidth=2, marker='o', s=20,
         transform=ccrs.PlateCarree(),
         )

plt.tight_layout()

plt.title('How to check if a point with a latitude and longitude is in a given area or not using shapely in python ?' ,fontsize=14)

plt.savefig("shapely_point_within_02.png", bbox_inches='tight', dpi=100)

plt.show()

plt.close()

However if we replace PlateCarree

poly = mpatches.Polygon(poly_corners, closed=True, fill=False, lw=1, fc=None, transform=ccrs.PlateCarree())

by Geodetic

poly = mpatches.Polygon(poly_corners, closed=True, fill=False, lw=1, fc=None, transform=ccrs.Geodetic())

we can see that the point is within the polygon !

Please take caution when selecting a coordinate system, as it can greatly impact your results.

From EPSG:4326 to EPSG:32618

from shapely.ops import transform

import pyproj

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

polygon = transform(project.transform, polygon)
point = transform(project.transform, point)

print(polygon.contains(point)) # check if polygon contains point
print(point.within(polygon)) # check if a point is in the polygon

returns

True
True

References