Python – Measuring Distance from Within Boundary Using Shapely

boundariespythonshapely

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:

exhibit a

exhibit b

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[0] for po in contracted for i in po.exterior.coords]
        int_y = [i[1] for po in contracted for i in po.exterior.coords]
    else:
        int_x = [i[0] for i in contracted.exterior.coords]
        int_y = [i[1] for i in contracted.exterior.coords]

    if type(dilated) is MultiPolygon:
        print("Uh oh. Dilated for " + polygon_file + " invalid.")
        ext_x = [i[0] for po in dilated for i in po.exterior.coords]
        ext_y = [i[1] for po in dilated for i in po.exterior.coords]
    else:
        ext_x = [i[0] for i in dilated.exterior.coords]
        ext_y = [i[1] 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

Best Answer

So myself and the PhD student I am working with figured out quite an elegant solution for what we were trying to do, which to remind you is that we wanted to acquire points that were within 500 of the boundary of the shape provided. Additionally, I also had to acquire those points NOT near the boundary that are within the shape - call these points as being within the 'core' of the shape.

What I did was this: 

def get_margin_core_points(all_points, shape, line, my_distance):
    margin_point_list, core_point_list = [], []

    for pt in all_points:
        x_avg, y_avg = np.mean([pt[0], pt[1]]), np.mean([pt[2], pt[3]])

        p = Point(x_avg, y_avg)

        # if within 500 of the boundary, on the margin.
        if p.distance(line) <= my_distance:
            margin_point_list.append(pt)
        # if that failed, try within the shape. Two comparisons only.
        elif p.within(shape):
            core_point_list.append(pt)

    return margin_point_list, core_point_list

# Don't worry about getting the shape and line, they're simple constructors,
# nothing fancy like I tried above. In other words, the original shape need
# NOT be valid for this to work.
def get_points_on_margin_and_core(all_points, annotations_file):
    all_class_a, all_class_b = all_points

    shape = load_shape(annotations_file)
    line = load_linestring(annotations_file)

    margin_class_a, core_class_a = get_margin_core_points(all_class_a, shape, line, 500)
    margin_class_b, core_class_b = get_margin_core_points(all_class_b, shape, line, 500)

    return np.asarray(margin_class_a), np.asarray(core_class_a),
           np.asarray(margin_class_b), np.asarray(core_class_b)

Essentially I treat the boundary as a linestring object AND a polygon. One can then efficiently determine whether the point lies within the original shape, whether it is close to the boundary, or neither.

Tah-dah! Quite pleased with this. Hope this helps others with similar issues.

Related Question