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
To check if a polygon contains multiple points I would use matplotlib
contains_points
, documented here: http://matplotlib.org/api/path_api.html#matplotlib.path.Path.contains_pointsThis does one big call using a numpy array, this is why it is efficient. Note that you can pass a radius which in fact inflates or delates the polygon, you can also transform (projections...) before doing the check.
I just did this to test if the points are actually inside the triangle:
Now when you use Google Maps and the polygon is mapped onto spheric coordinates, the triangle gets deformed, a thing to keep in mind.
Anyway, plotting your data with kml in Gookle Earth does show the point outside of the triangle as well?!
Same appearance as in the matplotlib image, Point 1 is slighlty outside of the triangle, when plotted in Euclidean 2D-coordinates. For geometric computations in geo-coordinates check QGIS Python Console or GDAL/OGR Tools. Or you would use the google maps api, just as in the example, that is linked on this page, where the topic 2D-geometries vs. geodesic geometries is coverd.
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 inmpl_toolkit
: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).
Seems to get the correct answer.