I'm new, although I am reasonably experienced in Python.
My question concerns Shapely, a library for spatial data. I am trying to assign some points as being within some distance d of the boundary of the shape. The main problem is that these shapes are in some cases very, very nasty indeed, with self-overlap and intersection galore - some examples illustrate what I am dealing with:
My initial approach was to try constructing buffered objects that were +/-d from the original shape using the
buffer() method. But this doesn't appear to work for me - the shapes are either overlapping, completely all over the place, or way too far in to be sensible.
I figured that another approach was to simply get the
distance() to the shape. But this does not work for those points within the shape -
distance() merely returns 0.
I am wondering if there is any way of getting the distance of a point to the boundary of a shape from within that shape, or more generally a way to do what I am trying to do. Any suggestions?
An example Python snippet of the approach that seems to work best is as follows:
# A function that reads in a likely invalid shape and tries to return # a 'valid' one. def construct_shape(poly, distance, delta, update): difference = 0 if type(poly.buffer(distance)) is MultiPolygon: difference = update(difference, delta) while type(poly.buffer(distance + difference)) is MultiPolygon: difference = update(difference, delta) new_poly = poly.buffer(distance + difference) return new_poly.buffer(-difference) else: return poly.buffer(distance) for (file_i, polygon_file) in enumerate(sorted(glob.glob("./line/*.annotations"))): # Do some processing with the files, which contain the points that define the shapes # tags is a list of points extracted from the current file # ... poly = Polygon(tags) # Iterative increasing approach # dilated is +500 from the shape dilated = construct_shape(poly, 500, 1, lambda x, y: x + y) # contracted is -500 from the shape, or in this case -1000 from dilated. contracted = construct_shape(dilated, -1000, 1, lambda x, y: x - y) if type(contracted) is MultiPolygon: print("Uh oh. Contracted for " + polygon_file + " invalid.") int_x = [i for po in contracted for i in po.exterior.coords] int_y = [i for po in contracted for i in po.exterior.coords] else: int_x = [i for i in contracted.exterior.coords] int_y = [i for i in contracted.exterior.coords] if type(dilated) is MultiPolygon: print("Uh oh. Dilated for " + polygon_file + " invalid.") ext_x = [i for po in dilated for i in po.exterior.coords] ext_y = [i for po in dilated for i in po.exterior.coords] else: ext_x = [i for i in dilated.exterior.coords] ext_y = [i for i in dilated.exterior.coords] # Plot the results plt.plot(ext_x, ext_y, color="r") plt.plot(int_x, int_y, color="b") patch = PolygonPatch(poly, fc="red", alpha=0.5, zorder=2) ax.add_patch(patch) plt.show()
The above is my current approach. What I had in mind for it would be something like this:
inner_margin = contracted outer_margin = dilated point = Point(...) if point.within(outer_margin) and not point.within(inner_margin): # do some stuff with point else: # ignore
For the other approach regarding points being some distance from the margin:
point = Point(...) # Assume we just got boundary from Polygon(tags) if point.distance(boundary) <= 500: # do some stuff with point else: # ignore