I am testing the point-in-polygon function with matplotlib and shapely.
Here is a map contains a Bermuda triangle polygon.
Google maps's point-in-polygon functions clearly shows testingPoint and testingPoint2 are inside of the polygon which is a correct result.
if I test the two points in matplotlib and shapely, only point2 passes the test.
In [1]: from matplotlib.path import Path
In [2]: p = Path([[25.774252, -80.190262], [18.466465, -66.118292], [32.321384, -64.75737]])
In [3]: p1=[27.254629577800088, -76.728515625]
In [4]: p2=[27.254629577800088, -74.928515625]
In [5]: p.contains_point(p1)
Out[5]: 0
In [6]: p.contains_point(p2)
Out[6]: 1
shapely shows the same result as matplotlib does.
In [1]: from shapely.geometry import Polygon, Point
In [2]: poly = Polygon(([25.774252, -80.190262], [18.466465, -66.118292], [32.321384, -64.75737]))
In [3]: p1=Point(27.254629577800088, -76.728515625)
In [4]: p2=Point(27.254629577800088, -74.928515625)
In [5]: poly.contains(p1)
Out[5]: False
In [6]: poly.contains(p2)
Out[6]: True
What is actually going on here? Is Google's algorithm better than those two?
Thanks
There are basically two ways of conducting Point in Polygon queries in Shapely: using a function called . within() that checks if a point is within a polygon. using a function called .
How to check if a point is inside a polygon in Python. To perform a Point in Polygon (PIP) query in Python, we can resort to the Shapely library's functions . within(), to check if a point is within a polygon, or . contains(), to check if a polygon contains a point.
Draw a horizontal line to the right of each point and extend it to infinity. Count the number of times the line intersects with polygon edges. A point is inside the polygon if either count of intersections is odd or point lies on an edge of polygon. If none of the conditions is true, then point lies outside.
Remember: the world isn't flat! If Google Maps' projection is the answer you want, you need to project the geographic coordinates onto spherical Mercator to get a different set of X and Y coordinates. Pyproj can help with this, just make sure you reverse your coordinate axes before (i.e.: X, Y or longitude, latitude).
import pyproj
from shapely.geometry import Polygon, Point
from shapely.ops import transform
from functools import partial
project = partial(
pyproj.transform,
pyproj.Proj(init='epsg:4326'),
pyproj.Proj('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs'))
poly = Polygon(([-80.190262, 25.774252], [-66.118292, 18.466465], [-64.75737, 32.321384]))
p1 = Point(-76.728515625, 27.254629577800088)
# Old answer, using long/lat coordinates
poly.contains(p1) # False
poly.distance(p1) # 0.01085626429747994 degrees
# Translate to spherical Mercator or Google projection
poly_g = transform(project, poly)
p1_g = transform(project, p1)
poly_g.contains(p1_g) # True
poly_g.distance(p1_g) # 0.0 meters
Seems to get the correct answer.
Although you have already accepted an answer, but in addition to @MikeT's answer I will add this for future visitors who might want to do the same with matplotlib
and basemap in mpl_toolkit
:
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
# Mercator Projection
# http://matplotlib.org/basemap/users/merc.html
m = Basemap(projection='merc', llcrnrlat=-80, urcrnrlat=80,
llcrnrlon=-180, urcrnrlon=180, lat_ts=20, resolution='c')
# Poly vertices
p = [[25.774252, -80.190262], [18.466465, -66.118292], [32.321384, -64.75737]]
# Projected vertices
p_projected = [m(x[1], x[0]) for x in p]
# Create the Path
p_path = Path(p_projected)
# Test points
p1 = [27.254629577800088, -76.728515625]
p2 = [27.254629577800088, -74.928515625]
# Test point projection
p1_projected = m(p1[1], p1[0])
p2_projected = m(p2[1], p2[0])
if __name__ == '__main__':
print(p_path.contains_point(p1_projected)) # Prints 1
print(p_path.contains_point(p2_projected)) # Prints 1
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