Shapely: Add functionality to drop 3rd dimension

[LostFan123]

There is an old question on GIS Stack Exchange about converting 3D geometries to 2D: Convert 3D WKT to 2D Shapely Geometry. I think this functionality should be included in the Shapely, so we could use it, for example, like:

>>> from shapely.geometry import Polygon
>>> p = Polygon([(0, 0, 0), (1, 0, 0), (1, 1, 0)])
>>> p.wkt
'POLYGON Z ((0 0 0, 1 0 0, 1 1 0, 0 0 0))'
>>> p2 = p.drop_z
>>> p2.wkt
'POLYGON ((0 0, 1 0, 1 1, 0 0))'

I've seen in one of the answers and in the docs that _this operation is not necessary_ as the 3rd dimension _has no effect on geometric analysis_. But when implementing my own function to determine left side of a split geometry (this is not yet implemented: https://github.com/Toblerity/Shapely/issues/589), I saw that this functionality could come handy:

from shapely.geometry import (LinearRing, 
                              LineString, 
                              Polygon)


def is_left(polygon: Polygon,
            line: LineString) -> bool:
    """
    Determines if the polygon is on the left side of the line
    according to:
    https://stackoverflow.com/questions/50393718/determine-the-left-and-right-side-of-a-split-shapely-geometry
    """
    ring = LinearRing([*line.coords, *polygon.centroid.coords])
    return ring.is_ccw

This code will fail for 3D geometries:

p = Polygon([(0, 0, 0), (1, 0, 0), (1, 1, 0)])
l = LineString([(0, 0, 0), (1, 0, 0)])
is_left(p, l)

will give this error:

AttributeError                            Traceback (most recent call last)
~/miniconda3/lib/python3.7/site-packages/shapely/speedups/_speedups.pyx in shapely.speedups._speedups.geos_linearring_from_py()

AttributeError: 'list' object has no attribute '__array_interface__'

During handling of the above exception, another exception occurred:

IndexError                                Traceback (most recent call last)
<ipython-input-55-555b9e2533fa> in <module>()
----> 1 is_left(p, l)

<ipython-input-52-7fad75b19ce3> in is_left(polygon, line)
      6     https://stackoverflow.com/questions/50393718/determine-the-left-and-right-side-of-a-split-shapely-geometry
      7     """
----> 8     ring = LinearRing([*line.coords, *polygon.centroid.coords])
      9     return ring.is_ccw

~/miniconda3/lib/python3.7/site-packages/shapely/geometry/polygon.py in __init__(self, coordinates)
     51         BaseGeometry.__init__(self)
     52         if coordinates is not None:
---> 53             self._set_coords(coordinates)
     54 
     55     @property

~/miniconda3/lib/python3.7/site-packages/shapely/geometry/polygon.py in _set_coords(self, coordinates)
     66     def _set_coords(self, coordinates):
     67         self.empty()
---> 68         ret = geos_linearring_from_py(coordinates)
     69         if ret is not None:
     70             self._geom, self._ndim = ret

~/miniconda3/lib/python3.7/site-packages/shapely/speedups/_speedups.pyx in shapely.speedups._speedups.geos_linearring_from_py()

IndexError: tuple index out of range

This is due to the fact that centroid of a Polygon is always returned in 2D (https://github.com/Toblerity/Shapely/issues/554), and a LinearRing can't be constructed from points having a different number of dimensions.

If there was a drop_z method, I would just write ring = LinearRing([*line.drop_z.coords, *polygon.centroid.coords]) instead of cluttering the code with things like line = LineString([xy[:2] for xy in list(line.coords)]) or implementing a function for that. Or even better, I would drop redundant 3rd dimension consisting only of zeros from the original parent polygon that I read from a file on the top level, so all the child geometries would have only 2 dimensions.

---

Shapely version: 1.6.4.post1, installed from conda.

[sgillies]

@LostFan123 I think a method to drop Z values would be an appropriate addition for 1.7. Thanks for suggesting it.

[astrojuanlu]

For other people arriving here from Google, I found a very simple method:

def _to_2d(x, y, z):
    return tuple(filter(None, [x, y]))

new_shape = shapely.ops.transform(_to_2d, shape)

(credit to @feenster and @hunt3ri from https://github.com/hotosm/tasking-manager/blob/master/server/services/grid/grid_service.py :heart:)

[sgillies]

I'm going to close this issue in favor of the 2-line solution above. Thanks @Juanlu001 !