shapely parallel_offset sometimes does not generate closed ring
Asked Answered
S

2

7

I am using the parallel_offset function of the shapely package to get offset structures to some polygons that are closed rings. I have several polygons at once, many with similar shapes. Around 10-25% of them, however, do not generate a closed ring from the parallel_offset. Here is a MWE of a shape that does not work:

import matplotlib.pyplot as plt
from shapely.geometry.polygon import LinearRing

def plot_line(ax, ob, color):
    x, y = ob.xy
    ax.plot(x, y, color=color, alpha=0.7, linewidth=3, 
            solid_capstyle='round', zorder=2)

polygon = [[-29.675, -30.675],
           [-28.4094, -29.4094],
           [-28.325, -29.325],
           [-28.325, -29.764],
           [-28.325, -29.7933],
           [-28.4587, -29.8274],
           [-28.4676, -29.8297],
           [-28.5956, -29.8814],
           [-28.6041, -29.8848],
           [-28.724, -29.953],
           [-28.732, -29.9576],
           [-28.8417, -30.0413],
           [-28.849, -30.0469],
           [-28.9466, -30.1445],
           [-28.9531, -30.151],
           [-29.0368, -30.2607],
           [-29.0424, -30.268],
           [-29.1106, -30.3879],
           [-29.1152, -30.3959],
           [-29.1669, -30.5239],
           [-29.1703, -30.5324],
           [-29.2044, -30.6661],
           [-29.2067, -30.675],
           [-29.6457, -30.675],
           [-29.675, -30.675]]

poly_line = LinearRing(polygon)
poly_line_offset = poly_line.parallel_offset(0.05, side="left", resolution=16, 
                                             join_style=2, mitre_limit=1)
fig = plt.figure()
ax = fig.add_subplot(111)
plot_line(ax, poly_line, "blue")
plot_line(ax, poly_line_offset, "green")
plt.show()

Blue plot original polygon, green plot offset polygon

As you can see, the green offset polygon does not close at the point that is first/last in the list of vertices. Other very similar shapes, however, do work as intended. They have the same data structure and also have the same start/end point, as does my example above. The join_style attribute does not change the outcome to what I want. Changing the resolution or distance does not help either. The documentation also does not address this issue.

Do you have any guidance? I am using shapely 1.6.3.

Sidonius answered 26/1, 2018 at 13:12 Comment(1)
The bug report has been filed hereHent
L
5

not completely sure why this happens, nevertheless you might use a workaround based on the buffer method:

poly_line = LinearRing(polygon)
poly_line_offset = poly_line.buffer(0.05, 
                       resolution=16, join_style=2, mitre_limit=1).exterior

With your data, this produces the (probably) desired result: enter image description here

Laurinelaurita answered 26/1, 2018 at 16:13 Comment(6)
Thanks for the answer! However, I found that buffer does not work in exactly the same way as parallel_offset does. I have a large set of (partially) nested contours. parallel_offset correctly recognizes what is "inside" and what is "outside", while buffer always at least misinterprets one contour.Sidonius
I also just noticed, despite using mitred join style, the bottom left corner is rounded. I guess this is a bug?Sidonius
And one last thing I noticed: buffer does not take direction as input like parallel_offset does. Negative values for distance do not work either.Sidonius
@Sidonius it seems to be something in the GEOSOffsetCurve method which shapely internally calls - I tried with PostGIS, and ST_OffsetCurve has the same behaviorLaurinelaurita
Where would I best report this issue?Sidonius
@ReblochonMasque I opened an issue there, we will see what happens next. I will report news here as edit to the OP.Sidonius
I
1

Here's a work around I did in my code. I basically rolled the LinearRing (shifting the start point along the ring), applied two offsets, and then added them back together. It's probably not an ideal solution, but hopefully can work as a starting point:

from shapely import ops, geometry
import numpy as np

# test geo:
ring_coords = [(0,0.1),(0,2),(4,2),(4,0)]
ring = geometry.LinearRing(ring_coords)


# shifts the ring by one point
rolled = LinearRing(np.roll(ring.coords[:-1], 2)) 

# apply the offsets
offset_ring = ring.parallel_offset(-0.2, side='right', resolution=3, join_style=2, mitre_limit=3)
offset_rolled = rolled.parallel_offset(-0.2, side='right', resolution=3, join_style=2, mitre_limit=3)

# combine the points
# assuming you started with two rings, backward should be empty
forward, backward = ops.shared_paths(offset_rolled, offset_ring) 
combined = geometry.LinearRing(ops.linemerge(forward))

Isfahan answered 19/10, 2021 at 2:0 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.