Check if geo-point is inside or outside of polygon
coordinate-systems, point-in-polygon, polygon, python, shapely
Solution
Here is a possible solution to my problem.
- Geographical coordinates must be stored properly. Example `np.array([[Lon_A, Lat_A], [Lon_B, Lat_B], [Lon_C, Lat_C]])`
- Create the polygon
- Create the point to be tested
- Use `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.
Special case: If point lies on borders of Polygon
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
Problem
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: - points are connected to each other by Great Circles! - I have researched in the internt and I ended up finding some similar questions like this one but I had no success since they all use of `.shp` files which I do not have.