4

I have the WKT of many linestrings which contains three dimensions, however my script only requires a two dimensional linestring and I therefore want to discard the third dimension.

My current code, using Shapely, does flatten from 3D to 2D, but it relies on the geometry's wkt method outputting as 2D (even if it is 3D). As this Shapely behaviour may change in the future, can anyone recommend a more robust method of flattening to 2D?

  • Python 2.7.13
  • Shapely 1.2.17
  • Ubuntu 12.10
>>> from shapely.geometry import LineString  
>>> from shapely.wkt import loads    
>>> line_3d_wkt = "LINESTRING (1 2 0, 4 3 0, 8 9 0)"   
>>> line_3d = loads(line_3d_wkt)
>>> line_3d.has_z
True
>>> list(line_3d.coords)
[(1.0, 2.0, 0.0), (4.0, 3.0, 0.0), (8.0, 9.0, 0.0)]  
>>> line_3d.wkt  # outputs as 2D even if 3D geometry   
'LINESTRING (1.0000000000000000 2.0000000000000000, 4.0000000000000000 3.0000000000000000, 8.0000000000000000 9.0000000000000000)'

My current technique to flatten the 3D geometry to 2D is:

>>> line_2d = loads(line_3d.wkt)
>>> line_2d.has_z
False
>>> list(line_2d.coords)
[(1.0, 2.0), (4.0, 3.0), (8.0, 9.0)]
3

With the new Shapely version of mwtoews following your link Set default WKB output dimensions to 3 for geometries with Z-dimension

>>> line_3d_wkt = "LINESTRING (1 2 0, 4 3 0, 8 9 0)"
>>> line_3d = loads(line_3d_wkt)
>>> line_3d.wkt
'LINESTRING Z (1 2 0, 4 3 0, 8 9 0)'
>>> line_2d = LineString([xy[0:2] for xy in list(line_3d.coords)]) 
>>> line_2d.wkt
'LINESTRING (1 2, 4 3, 8 9)'
2

I had the exact same problem as Rusty Magoo, but for multipolygons and polygons. Therefore I decided to build on answer of gene and wrote the following function, which works for all types of geometries.

from shapely.geometry import *

def remove_third_dimension(geom):
    if geom.is_empty:
        return geom

    if isinstance(geom, Polygon):
        exterior = geom.exterior
        new_exterior = remove_third_dimension(exterior)

        interiors = geom.interiors
        new_interiors = []
        for int in interiors:
            new_interiors.append(remove_third_dimension(int))

        return Polygon(new_exterior, new_interiors)

    elif isinstance(geom, LinearRing):
        return LinearRing([xy[0:2] for xy in list(geom.coords)])

    elif isinstance(geom, LineString):
        return LineString([xy[0:2] for xy in list(geom.coords)])

    elif isinstance(geom, Point):
        return Point([xy[0:2] for xy in list(geom.coords)])

    elif isinstance(geom, MultiPoint):
        points = list(geom.geoms)
        new_points = []
        for point in points:
            new_points.append(remove_third_dimension(point))

        return MultiPoint(new_points)

    elif isinstance(geom, MultiLineString):
        lines = list(geom.geoms)
        new_lines = []
        for line in lines:
            new_lines.append(remove_third_dimension(line))

        return MultiLineString(new_lines)

    elif isinstance(geom, MultiPolygon):
        pols = list(geom.geoms)

        new_pols = []
        for pol in pols:
            new_pols.append(remove_third_dimension(pol))

        return MultiPolygon(new_pols)

    elif isinstance(geom, GeometryCollection):
        geoms = list(geom.geoms)

        new_geoms = []
        for geom in geoms:
            new_geoms.append(remove_third_dimension(geom))

        return GeometryCollection(new_geoms)

    else:
        raise RuntimeError("Currently this type of geometry is not supported: {}".format(type(geom)))

Here are some unit-tests:

from unittest import TestCase

from shapely.geometry import *

from utils import remove_third_dimension


class TestGeometry(TestCase):
    def setUp(self):
        self.ext1 = LinearRing([(0, 0, 7), (5, 0, 7), (5, 5, 7), (0, 5, 7)])
        self.ext2 = LinearRing([(100, 100, 100), (101, 100, 100), (101, 101, 100), (100, 101, 100)])
        self.int1 = LinearRing([(3, 1, 3), (4, 1, 3), (4, 4, 3), (3, 4, 3)])
        self.int2 = LinearRing([(1, 1, 5), (2, 1, 5), (2, 2, 5), (1, 2, 5)])
        self.int3 = LinearRing([(1, 3, 5), (2, 3, 5), (2, 4, 5), (1, 4, 5)])

        self.line1 = LineString([(1, 1, 5), (2, 1, 5), (2, 2, 5), (1, 2, 5)])
        self.line2 = LineString([(8, 8, 5), (9, 8, 5), (9, 9, 5), (8, 9, 5)])
        self.empty_line = LineString()

        self.point1 = Point([0, 0, 100])
        self.point2 = Point([1, 1, 100])

    def test_polygon(self):
        pol1 = Polygon(self.ext1)
        self.assertAlmostEqual(pol1.area, 25, delta=1e-9)
        self.assertEqual(pol1.has_z, True)
        ext1 = pol1.exterior
        point1 = ext1.coords[0]
        self.assertEqual(len(point1), 3)
        pol1 = remove_third_dimension(pol1)
        self.assertAlmostEqual(pol1.area, 25, delta=1e-9)
        self.assertEqual(pol1.has_z, False)
        ext1 = pol1.exterior
        point1 = ext1.coords[0]
        self.assertEqual(len(point1), 2)

    def test_multipolygon(self):
        pol1 = Polygon(self.ext1)
        pol2 = Polygon(self.ext2)
        self.assertAlmostEqual(pol2.area, 1, delta=1e-9)
        self.assertEqual(pol2.has_z, True)
        multipol1 = MultiPolygon([pol1, pol2])
        self.assertAlmostEqual(multipol1.area, 26, delta=1e-9)
        self.assertEqual(multipol1.has_z, True)
        multipol1 = remove_third_dimension(multipol1)
        self.assertAlmostEqual(multipol1.area, 26, delta=1e-9)
        self.assertEqual(multipol1.has_z, False)

    def test_polygon_with_interior(self):
        pol3 = Polygon(self.ext1, [self.int1, self.int2, self.int3])
        self.assertAlmostEqual(pol3.area, 20, delta=1e-9)
        self.assertEqual(pol3.has_z, True)
        pol3 = remove_third_dimension(pol3)
        self.assertAlmostEqual(pol3.area, 20, delta=1e-9)
        self.assertEqual(pol3.has_z, False)

    def test_line(self):
        self.assertAlmostEqual(self.line1.length, 3, delta=1e-9)
        self.assertEqual(self.line1.has_z, True)
        line1 = remove_third_dimension(self.line1)
        self.assertAlmostEqual(line1.length, 3, delta=1e-9)
        self.assertEqual(line1.has_z, False)

    def test_empty_line(self):
        empty_line = remove_third_dimension(self.empty_line)
        self.assertEqual(empty_line.is_empty, True)

    def test_point(self):
        self.assertEqual(self.point1.has_z, True)
        point1 = remove_third_dimension(self.point1)
        self.assertEqual(point1.has_z, False)
        self.assertEqual(len(point1.coords[0]), 2)

    def test_multipoint(self):
        multipoint1 = MultiPoint([self.point1, self.point2])
        self.assertEqual(multipoint1.has_z, True)
        multipoint1 = remove_third_dimension(multipoint1)
        self.assertEqual(multipoint1.has_z, False)
        self.assertEqual(len(multipoint1.geoms), 2)

    def test_multiline(self):
        multiline1 = MultiLineString([self.line1, self.line2])
        self.assertEqual(multiline1.has_z, True)
        multiline1 = remove_third_dimension(multiline1)
        self.assertEqual(multiline1.has_z, False)
        self.assertEqual(len(multiline1.geoms), 2)
        self.assertAlmostEqual(multiline1.length, 6, delta=1e-9)

    def test_collection(self):
        col = GeometryCollection([self.ext2, self.line2, self.point2])
        self.assertEqual(col.has_z, True)
        col = remove_third_dimension(col)
        self.assertEqual(col.has_z, False)
        self.assertEqual(len(col.geoms), 3)
1

I do not understand your question, because you offer your solution.

In addition, this operation is not necessary, since:

A third z coordinate value may be used when constructing instances, but has no effect on geometric analysis. All operations are performed in the x-y plane.

from The Shapely User Manual

>>> line_3d.has_z
True
>>> buf = line_3d.buffer(5)
>>> buf.has_z
False
  • Primarily I a want to remove the third dimension as my application will be storing many thousands of geometry objects for later processing and I am presuming that the removal of this third dimension will be beneficial in terms of performance and memory consumption. It appears from the Shapely gitgub site link, that the technique I am currently using may well not be available in the future. – Rusty Magoo Jul 28 '13 at 10:06
  • 2
    Then, use list slicing: p= LineString([x[0:2] for x in list(line_3d.coords)]) – gene Jul 28 '13 at 10:25
  • This doesn't work anymore :( – astrojuanlu Aug 23 at 16:48
1

The functionality to drop 3rd dimension will be probably implemented in Shapely 1.7. See related GitHub issue.
While it's not ready yet, I'm using the following function that supports all the geometries (Python >= 3.7). It is much more concise than the one from Rauni's answer (passes all their unit-tests as well).

from functools import singledispatch
from operator import itemgetter
from typing import TypeVar

from shapely.geometry import Polygon
from shapely.geometry.base import (BaseGeometry,
                                   BaseMultipartGeometry)

Geometry = TypeVar('Geometry', bound=BaseGeometry)


@singledispatch
def drop_z(geometry: Geometry) -> Geometry:
    """
    Removes Z-coordinate from a geometry object.
    Won't be necessary after it will be implemented in Shapely:
    https://github.com/Toblerity/Shapely/issues/709
    """
    if geometry.is_empty:
        return geometry
    geometry_type = type(geometry)
    xy_coordinates = map(itemgetter(0, 1), geometry.coords)
    return geometry_type(list(xy_coordinates))


@drop_z.register
def _(geometry: BaseMultipartGeometry):
    geometry_type = type(geometry)
    return geometry_type(list(map(drop_z, geometry.geoms)))


@drop_z.register
def _(geometry: Polygon):
    if geometry.is_empty:
        return geometry
    new_exterior = drop_z(geometry.exterior)
    new_interiors = list(map(drop_z, geometry.interiors))
    return Polygon(new_exterior, new_interiors)

Simple example usage:

>>> print(drop_z(Point(1, 2, 3)))
POINT (1 2)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.