cga.py 1.7 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950
  1. """Shapely CGA algorithms."""
  2. import numpy as np
  3. import shapely
  4. def signed_area(ring):
  5. """Return the signed area enclosed by a ring in linear time.
  6. Algorithm used: https://web.archive.org/web/20080209143651/http://cgafaq.info:80/wiki/Polygon_Area
  7. """
  8. coords = np.array(ring.coords)[:, :2]
  9. xs, ys = np.vstack([coords, coords[1]]).T
  10. return np.sum(xs[1:-1] * (ys[2:] - ys[:-2])) / 2.0
  11. def _reverse_conditioned(rings, condition):
  12. """Return a copy of the rings potentially reversed depending on `condition`."""
  13. condition = np.asarray(condition)
  14. if np.all(condition):
  15. rings = shapely.reverse(rings)
  16. elif np.any(condition):
  17. rings = np.array(rings)
  18. rings[condition] = shapely.reverse(rings[condition])
  19. return rings
  20. def _orient_polygon(geometry, exterior_cw=False):
  21. if geometry is None:
  22. return None
  23. if geometry.geom_type in ["MultiPolygon", "GeometryCollection"]:
  24. return geometry.__class__(
  25. [_orient_polygon(geom, exterior_cw) for geom in geometry.geoms]
  26. )
  27. # elif geometry.geom_type in ["LinearRing"]:
  28. # return reverse_conditioned(geometry, is_ccw(geometry) != ccw)
  29. elif geometry.geom_type == "Polygon":
  30. rings = np.array([geometry.exterior, *geometry.interiors])
  31. reverse_condition = shapely.is_ccw(rings)
  32. reverse_condition[0] = not reverse_condition[0]
  33. if exterior_cw:
  34. reverse_condition = np.logical_not(reverse_condition)
  35. if np.any(reverse_condition):
  36. rings = _reverse_conditioned(rings, reverse_condition)
  37. return geometry.__class__(rings[0], rings[1:])
  38. return geometry
  39. _orient_polygons_vectorized = np.frompyfunc(_orient_polygon, nin=2, nout=1)