cuspatial icon indicating copy to clipboard operation
cuspatial copied to clipboard

[FEA] Integrate `GeoSeries` with `read_polygon_shapefile`

Open thomcom opened this issue 2 years ago • 3 comments

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.

thomcom avatar Jul 29 '22 16:07 thomcom

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

zhangjianting avatar Aug 03 '22 23:08 zhangjianting

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:

  1. Simply declaring not supporting multi-polygons or GeometryCollection.
  2. Modify the Shapefile reader and extend the code of all relevant operations to truly support multi-polygons.

zhangjianting avatar Aug 04 '22 00:08 zhangjianting

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.

zhangjianting avatar Aug 25 '22 15:08 zhangjianting

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.

thomcom avatar Aug 25 '22 19:08 thomcom

We should definitely have an area computation as well, and it looks quite simple. I'll file an issue for it.

thomcom avatar Aug 25 '22 19:08 thomcom

Ready for your re-review @harrism @isVoid

thomcom avatar Aug 29 '22 21:08 thomcom

@harrism @isVoid this is cleaner and readier for another round of review.

thomcom avatar Sep 07 '22 23:09 thomcom

@gpucibot merge

thomcom avatar Sep 19 '22 18:09 thomcom