[GIS] How to convert a geometry into LineString whose start and end point is equal

geometryshapely

I'm using shapely to calculate the boundary of individual LineString from a MultiLineString. The code was working fine until I encountered a case where the extracted LineString has the same start and end point. Though it should be a Polygon but upon printing the wkt version of the shapely object I found out that shapely is considering it as a LineString. Now the boundary method is failing here as I'm getting MULTIPOINT EMPTY as output. Is this a bug in shapely or I'm doing something wrong ?

Here's my code:

for lines in multilines:
    print lines.wkt
    boundary_points = lines.boundary

The output is:

LINESTRING (6.799200988255513 51.44578708679462, 6.79938 51.4457, 
6.80022 51.44523, 6.80187 51.44439, 6.80447 51.44315, 6.80474 51.44293, 6.80487 51.44278, 6.80494 51.44264, 6.80496 51.44246, 
6.80481 51.44219, 6.80467 51.44205, 6.8045 51.44196, 6.80427 51.44188, 6.80408 51.44184, 6.80387 51.44182, 6.80363 51.44183, 
6.80327 51.44187, 6.80294 51.44195, 6.80256 51.4421, 6.80191 51.44245, 6.80136 51.44287, 6.80095 51.44323, 6.80017 51.44407, 
6.79988 51.44447, 6.79963 51.44488, 6.7994 51.44532, 6.79921 51.44575, 6.799200988255513 51.44578708679462)

Best Answer

A LineString with the same start and end point is not a Polygon but a LinearRing. You need to understand the difference. If you want a Polygon from the LinearRing

# test if your output is a LinearRing
line.is_ring
True
# thus you can make a Polygon with the LinearRing
poly = Polygon(yourline)

Now from the manual of Shapely

A ring’s boundary is empty.

But poly.boundary not why ? Because the poly.boundary = the LinearRing of the Polygon

ring = LinearRing([(0, 0), (1, 1), (1, 0)])
for pt in list(ring.coords):
    print Point(pt)
 POINT (0 0)
 POINT (1 1)
 POINT (1 0)
 POINT (0 0)
poly = Polygon(ring)
print ring.boundary
MULTIPOINT EMPTY
print poly.boundary # the result is a LineString and not Points
LINESTRING (0 0, 1 1, 1 0, 0 0)
for i in list(poly.exterior.coords):
   print Point(i)
POINT (0 0)
POINT (1 1)
POINT (1 0)
POINT (0 0)

Related Solutions

Python Shapely – How to Sort Points in a MultiLineString for Continuous LineString

You might be able to just use shapely.ops.linemerge. I haven't used it much, and cannot speak to how it handles potential floating point differences where the geometries in the MultiLineString meet.

from shapely import wkt
from shapely.ops import linemerge


mls = wkt.loads('MULTILINESTRING ((120 105, 140 98),(160 103, 140 98),(100 100, 120 105))')

line_string = linemerge(mls.geoms)
# LINESTRING (100 100, 120 105, 140 98, 160 103)

Python OGR to Shapely – How to Convert Geometries with M-Values to Those Without

Here's a potential solution following suggestions from jbalk's and user30184's comments. It uses OGR's SetMeasured function to drop the M-values and relies on WKB (instead of WKT) to do the actual conversion between OGR and Shapely.

def from_ogr_to_shapely(ogr_geom):
    # Creating a copy of the input OGR geometry. This is done in order to 
    # ensure that when we drop the M-values, we are only doing so in a 
    # local copy of the geometry, not in the actual original geometry. 
    ogr_geom_copy = ogr.CreateGeometryFromWkb(ogr_geom.ExportToIsoWkb())
    
    # Dropping the M-values
    ogr_geom_copy.SetMeasured(False)

    # Generating a new shapely geometry
    # Note: I needed to typecast the results from
    # `ogr_geom_copy.ExportToIsoWkb()` to bytes since, in newer versions 
    # of GDAL/OGR and shapely, the outputs of that method were incompatible 
    # with the expected inputs for the `shapely.wkb.loads` method, 
    # leading to an error.
    shapely_geom = shapely.wkb.loads(bytes(ogr_geom_copy.ExportToIsoWkb()))

    return shapely_geom

As mentioned in the commented section in the code above, we need to create a local copy of the input OGR geometry. This is because the function used to drop the M-values (the SetMeasured function) actually modifies the original input geometry. So in order to ensure that our original geometry remains intact, we need to create a local copy. Then we can modify the local copy to our heart's content without fear of changes creeping upstream to the original OGR geometry.

Here's the function being applied to the inputs in the question:

print(from_ogr_to_shapely(pt_2D).wkt)
# POINT (0 0) # right!

print(from_ogr_to_shapely(line1_2D).wkt)
# LINESTRING (0 0, 1 1) # right!

print(from_ogr_to_shapely(mline_2D).wkt)
# MULTILINESTRING ((0 0, 1 1), (3 1, 4 1)) # right!

print(from_ogr_to_shapely(pt_3D).wkt)
# POINT Z (0 0 0) # right!

print(from_ogr_to_shapely(line1_3D).wkt)
# LINESTRING Z (0 0 0, 1 1 1) # right!

print(from_ogr_to_shapely(mline_3D).wkt)
# MULTILINESTRING Z ((0 0 0, 1 1 1), (3 1 3, 4 1 6)) # right!

Edit from the future

As pointed out by @ermaure in their answer, it seems like more recent versions of GDAL/OGR and shapely have changed things a bit and now the outputs of the ogr_geom.ExportToIsoWkb method are incompatible with the expected inputs for the shapely.wkb.loads method, causing the shapely.wkb.loads(ogr_geom_copy.ExportToIsoWkb()) statement to yield an error. This can be fixed by forcing the result of the ExportToIsoWkb method to be typecast to bytes. Explicitly, the command needs to change to shapely.wkb.loads(bytes(ogr_geom_copy.ExportToIsoWkb())). I have edited the code above to reflect that change.

Special thanks to @ermaure. Please be sure to upvote their answer as well!

Best Answer

Related Solutions

Python Shapely – How to Sort Points in a MultiLineString for Continuous LineString

Python OGR to Shapely – How to Convert Geometries with M-Values to Those Without

Edit from the future

Related Question