Check if geo-point is inside or outside of polygon

Question

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.

Oxyacetylene

Try converting this algorithm to Python wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html#The C Code – Heavierthanair

python does not have packages that do anything. it has a small number of pre built modules. packages are usually supplied by community. – Ardeb

Is the polygon always convex? – Hereon

In general no, it could also be concave – Oxyacetylene

Take a look at this question and implement this algorithm. Its not difficult #218078 – Hereon

Just in case: you can always cast a ray from your point to a middle point of any of the polygon's sides. If your ray crosses polygon's sides an even number of times, the point is on the outside. Works with convex and concave polygons; works on a sphere surface (and likely any 1-connected surface) using a geodesic for the ray. Has an edge case when a ray passes exactly through a vertex: you need to check whether the edges incident to the vertex are on the same side of the ray. – Erkan

Answer 1

Here is a possible solution to my problem.

1. Geographical coordinates must be stored properly. Example np.array([[Lon_A, Lat_A], [Lon_B, Lat_B], [Lon_C, Lat_C]])
2. Create the polygon
3. Create the point to be tested
4. 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 

Oxyacetylene

It is better to write the first latitude followed by longitude. Nothing wrong with the logic here, but to be safe than sorry. – Sori

@FedericoGentile - Shapely does not use great circle distances. It uses Euclidean – Eiland

Does this work for a region bounded by smooth curves as well? Not just a polygon. – Header

Thanks! The answer is compact but answer all the questions I need! Cool! – Deviation

Shapely also fails when the polygon intersects with the dateline – Wardrobe

Answer 2

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)

Tennyson

I've been having a lot of trouble installing shapely. This is a much cleaner solution - thank you! – Innuendo

Although boolean_point_in_polygon(point, polygon) perfectly addresses OP’s use case, points_within_polygon(points, polygon) is worth checking for aggregation use cases. – Honebein

Answer 3

Another way to do it is by using the even-odd algorithm explained in this link https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html The python code is given in wikipedia https://en.wikipedia.org/wiki/Even–odd_rule

Folks, just remember that the ORDER OF POINTS that make the polygon MATTER! I mean, different order results in different polygons.

Armorer

Answer 4

You can use pygeodesy package, it has no dependencies on system level geo packages, and it uses Kenneth Gade's n-vector approach,

https://github.com/mrJean1/PyGeodesy

Just pip install pygeodesy

Sample code

from pygeodesy.sphericalNvector import LatLon

p = LatLon(45.1, 1.1)
b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
print (p.isenclosedBy(b))

This should output True

Ireneirenic