cuspatial
cuspatial copied to clipboard
[FEA] Integrate `GeoSeries` with `read_polygon_shapefile`
This useful PR adds a GeoSeries
constructor for the tuple returned by read_polygon_shapefile
. Now users can say
geoseries = cuspatial.GeoSeries(cuspatial.read_polygon_shapefile('the_shapefile'))
and load the data directly into GeoArrow.
This PR touches many files because I add a reversed
argument to python and C++ so that the user can swap back and forth if they desire. This was required because the original implementation was reverse-ordering the polygons, which does not match GeoPandas results when reading a shapefile. Now it does, and is tested in both configurations.
There are very informative discussions on the ordering of polygon vertices: https://gis.stackexchange.com/questions/119150/order-of-polygon-vertices-in-general-gis-clockwise-or-counterclockwise The OGC SFS and ESRI shapefile have different orderings. As long as the inner rings and the outer rings have different directions, PIP algorithms should still work. Also see https://wrfranklin.org/Research/Short_Notes/pnpoly.html#Listing%20the%20Vertices
The Shapefile Reader API does have a problem in dealing with multi-polygons. The initial design was to break multi-polygons (as well as GeometryCollection) into multiple (single/regular) polygons (with or without holes). Yet the implementations was to consolidate the rings of multi-polygons, which essentially treats the first ring as the outer ring and the rest as the inner rings, i.e., a single/regular polygon with holes. Subsequently, it will break PIP algorithms. There are two solutions:
- Simply declaring not supporting multi-polygons or GeometryCollection.
- Modify the Shapefile reader and extend the code of all relevant operations to truly support multi-polygons.
There is a simple way to determine the orientation of a ring/curve, e.g., signed area. https://en.wikipedia.org/wiki/Curve_orientation Also https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order It actually might be useful to return the area of rings for "advanced" users or applications that require computing areas of polygons. A Thrust segmented reduction (prefix-sum) will do.
Ah nice! Computing the handedness of the polygon is nothing more than the determinant of the matrix formed by the first three points. That handedness is guaranteed to continue through all of the rest of the points. It is tempting to implement this for completeness.
We should definitely have an area computation as well, and it looks quite simple. I'll file an issue for it.
Ready for your re-review @harrism @isVoid
@harrism @isVoid this is cleaner and readier for another round of review.
@gpucibot merge