Chung - 4 months ago 25

Python Question

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

`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

Answer

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.

Source (Stackoverflow)

Comments