Comment déterminer si un point, spécifié par sa latitude et sa longitude, se situe dans une zone donnée en utilisant la bibliothèque shapely en Python ?

Utilisant la bibliothèque Shapely en Python, il est possible de vérifier si un point avec une latitude et une longitude se trouve dans une zone ou non. Tout d'abord, vous devez créer un objet de géométrie Shapely pour votre zone d'intérêt (polygone). Ensuite, vous pouvez créer un objet Point Shapely avec votre latitude et votre longitude. Enfin, vous pouvez vérifier si ce point se trouve dans la zone définie par le polygone en utilisant les méthodes contains ou within, exemples :

Installer Shapely

Shapely peut être installé, par exemple, avec geopandas:

pip install geopandas

Vérifiez si un point se trouve à l'intérieur d'un polygone donné.

Veuillez importer les bibliothèques Python requises :

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

import numpy as np

Définissez un point géométrique

Définissez un point cible, par exemple.

target_lat = 43.7580
target_lon = -65.3208

point = Point(target_lon,target_lat) # create point

Définissez un polygone

Définissez une zone en utilisant un polygone shapely.

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

Utilisation de la méthode contains() de shapely.

Maintenant, vous pouvez utiliser cette méthode pour vérifier si un point avec une latitude et une longitude se trouve dans une zone d'intérêt. Il est important de noter que les coordonnées doivent être fournies dans le même système de coordonnées que celui utilisé pour votre polygone, sinon cela produira des résultats erronés.

polygon.contains(point)

donne ici

True

Utiliser la méthode within() de shapely

Pour vérifier si un point est dans le polygone

point.within(polygon)

retourne ici

True

Visualisation des données.

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

Transformer le polygone Shapely et le point dans un autre système de coordonnées.

Maintenant, considérons l'exemple suivant :

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

va donner

False
False

Traçons les données.

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

Cependant, si nous remplaçons PlateCarree

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

par Geodetic

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

Nous pouvons voir que le point est à l'intérieur du polygone !

Veuillez faire preuve de prudence lors de la sélection d'un système de coordonnées, car cela peut avoir un impact considérable sur vos résultats.

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

donne

True
True

Références