3


What I’m 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, then can’t be used with concave-hull because there are no longer any points.

How can I make this work?

Looking for Concave_by_Attribute

  • For the Concave Hull step, there is an option to "Split multipart geometries into Single part". Did you enable this? – klewis Jul 10 '17 at 19:33
  • Thanks for the response Klewis. Yes I have tried that but on a multipart layer I get the error: float division by zero See log for more details – Phil_in_Tx Jul 10 '17 at 21:39
  • I'm getting the same error, QGIS 2.16.1 – klewis Jul 10 '17 at 22:53
  • it looks like there is a single light blue point in the dataset, that might be throwing the division by zero error, since it can't forma polygon with less than 3 points – TDavis Jul 11 '17 at 0:45
  • Excellent point TDavis. I thought you had it but after removing any attribute with less than 4 points I still get the same error. – Phil_in_Tx Jul 11 '17 at 1:17
2

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 :)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.