4

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?

Looking for Concave_by_Attribute

5
  • 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
3

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

1

You can use two options.

I take in my example a point layer named points (not multipart but will not change anything) with colum fid, mygroup (for the group info).

Option 1:

I add a Virtual layer with following SQL query

SELECT mygroup, convex_hull(st_union(geometry)) as geom, sum(num_geometries(geometry)) as tot
FROM points
GROUP BY mygroup
HAVING sum(num_geometries(geometry)) > 2

For human explanation, the query says group all geometry using the content from mygroup column, union them. From the union convert them to convex hulls. Only do it when the number of points in the group is at least 3 points with HAVING sum(num_geometries(geometry)) > 2

See the screen captures

Virtual layer windows

Main QGIS windows overview

Option 2:

You can use processing algorithm "Execute SQL" (in section "Vector General") as seen below

Overview of "Execute SQL" in Processing Toolbox

Then, fill the info like below to get an output similar to the one with virtual layer but in a new layer, temporary or a particular file depending of your output choice.

Execute SQL Windows

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.