omanges / turfpy

A Python library for performing geospatial data analysis which reimplements turf.js.
MIT License
121 stars 27 forks source link

Trouble with boolean_point_in_polygon() on irregular polygons. #100

Open mcgeochd opened 3 years ago

mcgeochd commented 3 years ago

I have a program that uses scipy spatial SpericalVoronoi to generate irregular non-overlapping polygons covering the surface of a sphere and I'm trying to use turfy to rasterize the polygons with boolean_point_in_polygon(). Some example outputs are:

https://imgur.com/a/nIQTSNF https://imgur.com/a/8DpitYe

The top image is a method of rasterizing that calculates the distance from a raster square to every edge of every polygon (an edge in this case refers to a single great circle segment between vertices) by projecting the point onto the line segment, finds the edge the raster point is closest to, and determines which side of this edge its on using the sign of the cross track distance (using a custom implementation, not the turpy implementation), and therefore which polygon it belongs to. This has some issues, as can be seen with the inlcusions of the wrong coloured polygon. This occurs when two edges are equally close to a raster point, i.e. the vertex shared by two edges is the closest point, and the wrong edge is selected, resulting in an incorrect rasterization.

I'm trying to get around this by using turfpy to rasterize the polygons, the results of which are shown in the lower image. The trouble is that the geojson polygon needs the polygon vertices to be supplied in clockwise order. While I can create a directed edge for each polygon, I'm having trouble detecting whether a given directed edge is defined clockwise or not, if it isn't I can just reverse the list order. I've tried using Python implementations of two methods described here without much luck: https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order/1180256#1180256

Signed Area

area = 0
for i in range(len(verts)):
    v1 = verts[i]
    v2 = verts[(i+1)%len(verts)]
    area += (v2[0] - v1[0]) * (v2[1] + v1[1])
if area < 0: verts = verts[::-1]
self.poly = Polygon([verts])

Convex Hull Point

minLat = 1e18
A = None
for v in verts:
    if v[1] < minLat:
        minLat = v[1]
        A = v
    elif v[1] == minLat and v[0] > A[0]:
        A = v
Aindex = indexOf(verts, A)
A = np.array(A)
B = np.array(verts[(Aindex-1)%len(verts)])
C = np.array(verts[(Aindex+1)%len(verts)])
AB = B - A
AC = C - A
cross = np.cross(AB, AC)
if cross < 0: verts = verts[::-1]
self.poly = Polygon([verts])

If anybody knows what is going wrong, or what algorithm turfpy uses for boolean_point_in_polygon so I can write an implementation that is direction agnostic, that would be super helpful.

mcgeochd commented 3 years ago

I fixed a problem where the GeoJson polygons have to have the same first and last point, which I had missed previously, but the output is still full of errors: https://imgur.com/a/lV2ZqsY