I am using python and I have defined the latitudes and longitudes (in degrees) of a polygon on the map. My goal is to check if a generic point P
of coordinates x,y
falls within such polygon. I would like therefore to have a function that allows me to check such condition and return True
or False
if the point is inside or outside the polygon.
In this example the point is outside so the result would be False
Question: Is there a library/package that allows to reach my goal? if yes which one do you recommend? would you be able to give a small example on how to use it?
Here is the code I have written so far:
import numpy as np # Define vertices of polygon (lat/lon) v0 = [7.5, -2.5] v1 = [2, 3.5] v2 = [-2, 4] v3 = [-5.5, -4] v4 = [0, -10] lats_vect = np.array([v0[0],v1[0],v2[0],v3[0],v4[0]]) lons_vect = np.array([v0[1],v1[1],v2[1],v3[1],v4[1]]) # Point of interest P x, y = -6, 5 # x = Lat, y = Lon ## START MODIFYING FROM HERE; DO NOT MODIFY POLYGON VERTICES AND DATA TYPE # Check if point of interest falls within polygon boundaries # If yes, return True # If no, return False
In order to plot the polygon and the point of interest I used cartopy and I wrote the following lines of code:
import cartopy.crs as ccrs import matplotlib.pyplot as plt ax = plt.axes(projection=ccrs.PlateCarree()) ax.stock_img() # Append first vertex to end of vector to close polygon when plotting lats_vect = np.append(lats_vect, lats_vect[0]) lons_vect = np.append(lons_vect, lons_vect[0]) plt.plot([lons_vect[0:-1], lons_vect[1:]], [lats_vect[0:-1], lats_vect[1:]], color='black', linewidth=1, transform=ccrs.Geodetic(), ) plt.plot(y, x, '*', # marker shape color='blue', # marker colour markersize=8 # marker size ) plt.show()
Note:
.shp
files which I do not have.Use polygon. contains(point) to test if point is inside ( True ) or outside ( False ) the polygon.
Algorithm: For a convex polygon, if the sides of the polygon can be considered as a path from any one of the vertex. Then, a query point is said to be inside the polygon if it lies on the same side of all the line segments making up the path.
A convex polygon is a simple polygon (with no self intersections) such that any line segment between two points inside of the polygon lies completely inside of it. The point will be inside a convex polygon if and only if it lies on the same side of the support line of each of the segments.
Here is a possible solution to my problem.
np.array([[Lon_A, Lat_A], [Lon_B, Lat_B], [Lon_C, Lat_C]])
polygon.contains(point)
to test if point is inside (True
) or outside (False
) the polygon.Here is the missing part of the code:
from shapely.geometry import Point from shapely.geometry.polygon import Polygon lons_lats_vect = np.column_stack((lons_vect, lats_vect)) # Reshape coordinates polygon = Polygon(lons_lats_vect) # create polygon point = Point(y,x) # create point print(polygon.contains(point)) # check if polygon contains point print(point.within(polygon)) # check if a point is in the polygon
Note: the polygon does not take into account great circles, therefore it is necessary to split the edges into many segments thus increasing the number of vertices.
E.g. print(Polygon([(0, 0), (1, 0), (1, 1)]).contains(Point(0, 0)))
will fail
So one can use
print(polygon.touches(point)) # check if point lies on border of polygon
There is also an emerging python library turfpy. which is used for geospatial analysis.
PyPI
Github
Example:
from turfpy.measurement import boolean_point_in_polygon from geojson import Point, Polygon, Feature point = Feature(geometry=Point((-46.6318, -23.5523))) polygon = Polygon( [ [ (-46.653, -23.543), (-46.634, -23.5346), (-46.613, -23.543), (-46.614, -23.559), (-46.631, -23.567), (-46.653, -23.560), (-46.653, -23.543), ] ] ) boolean_point_in_polygon(point, polygon)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With