[GIS] Concave Hull by Attribute type solution

Question

What I am looking for is a concave_hull-by-attribute type solution, or something that results in the same effect.
I have a multipart point layer that I want to group by attribute with minimum overlap (fig1). The Concave hull plugin does not work on multipart layers. It wants to group ALL the points on a layer together (fig2). Convex plug does accept hull-by-attribute but does not closely follow the points to the degree I need (fig3).

I have also tried breaking the multipart layer into separate layers using the `Spit` vector layer plugin but that process combines all the points to one in the attribute table. So, it cannot be used with concave-hull because there are no longer any points.

How can I make this work?

Answer 1

I developed following python-solution for my case. It should be easily applicable to your case, even without much python knowledge.

Import statements:

import numpy as np
import shapely
import fiona
import shapely.geometry as geometry
from shapely.ops import cascaded_union, polygonize
from scipy.spatial import Delaunay
from shapely.geometry import mapping
import math

Importing the points, where "feature_column_name" could be something like "class_id":

def load_points(shapefile, feature_column_name):
    shapefile = fiona.open(shapefile)
    points = np.array(
        [[p["geometry"]["coordinates"][0], p["geometry"]["coordinates"][1], p["properties"][feature_column_name]]
         for p in shapefile])
    return points

Then we use a modified version of the alpha-shape function from this tutorial. It now takes a buffer value as well. I have changed the meaning of the alpha, as I find it easier to use. Large values mean that your polygon contains more points. Also, we check if the resulting polygon contains at least 90% of the points. If not or if the number of points is below 20, we return the convex hull!

def alpha_shape(points, alpha, buffer=0):

    if len(points) < 20:
        return geometry.MultiPoint(points).convex_hull, None

    def add_edge(edges, edge_points, coords, i, j):
        """
        Add a line between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            return
        edges.add( (i, j) )
        edge_points.append(coords[ [i, j] ])

    coords = np.array(points)

    try:
        tri = Delaunay(coords)
    except:
        return None, None

    edges = set()
    edge_points = []
    # loop over triangles:
    # ia, ib, ic = indices of corner points of the
    # triangle
    for ia, ib, ic in tri.vertices:
        pa = coords[ia]
        pb = coords[ib]
        pc = coords[ic]
        # Lengths of sides of triangle
        a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
        b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
        c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
        # Semiperimeter of triangle
        s = (a + b + c)/2.0
        try:
            # Area of triangle by Heron's formula
            area = math.sqrt(s*(s-a)*(s-b)*(s-c))
        except ValueError:
            area = 0
        # print(ia, ib, ic, area)
        circum_r = a*b*c/(4.0*area+10**-5)
        # print(circum_r)
        # Here's the radius filter.
        if circum_r < alpha:
            add_edge(edges, edge_points, coords, ia, ib)
            add_edge(edges, edge_points, coords, ib, ic)
            add_edge(edges, edge_points, coords, ic, ia)
    try:
        m = geometry.MultiLineString(edge_points)
        triangles = list(polygonize(m))
        concave_hull = cascaded_union(triangles)
        concave_hull = concave_hull.buffer(buffer)

    except:
        return None, None

    # Lets check, if the resulting polygon contains at least 90% of the points.
    # If not, we return the convex hull.

    points_total = len(points)
    points_inside = 0

    for p in shapely.geometry.MultiPoint(points):
        points_inside += concave_hull.contains(p)

    if points_inside/points_total<0.9:
        return geometry.MultiPoint(points).convex_hull, None
    elif not concave_hull.is_empty:
        return concave_hull, edge_points
    else:
        return None, None

We can then loop over all our different points using the following function. "features" has to be iterable. You have to know the classes of your points in advance. So "features" could for example contain the class_ids [1,3,5,...] or just be range(1234).

def produce_alpha_polygons(features, alpha, buffer=0):
    alpha_polygons = []

    for i, f in enumerate(features):
        #print("{}/{}".format(i+1, len(features)))
        #select certain features
        current_points = points[points[:,-1]==f][:,:2]
        alpha_poly, _ = alpha_shape(current_points, alpha, buffer)
        alpha_polygons += [alpha_poly]

    return alpha_polygons

As we now have the polygons, we can export them:

def write_polygons(output_file, polygons):
    # Define a polygon feature geometry with one attribute
    schema = {
        'geometry': 'Polygon',
        'properties': {'id': 'int'},
    }

    # Write a new Shapefile
    with fiona.open(output_file, 'w', 'ESRI Shapefile',
                    schema) as f:
        for i, poly in enumerate(polygons):
            if not (isinstance(poly, shapely.geometry.Polygon) or isinstance(poly, shapely.geometry.MultiPolygon)): continue
            f.write({
                'geometry': mapping(poly),
                'properties': {'id': i},
            })

Usage:

points = load_points(input_shapefile, 'your attrib field name')
alpha_polygons = produce_alpha_polygons(range(1234), alpha=50, buffer=5)
write_polygons(output_file_name, alpha_polygons)

I used this method to produce concave hulls for roughly 2 million points in 6400 different clusters. I hope, this helps, even if I'm a year late :)