Concave Hulls ... the elusive container

14224
18
03-10-2018 04:14 PM
Labels (1)
DanPatterson_Retired
MVP Emeritus
7 18 14.2K

Concave hulls...

Tons of examples, some suited to some point clouds, others not so much.

It has greatly amused me over the years that people spend so much time trying to find the 'corner cases' where a particular implementation fails.  For example, the ever popular C - shaped object. 

I know... in optical fields, being able to identify the letters of the alphabet or curly snaked-shaped patterns is important.... but, how many times have you sent your field assistant to collect data points with weird shapes.

Toolbox links

Concave Hulls  this is a separate toolbox

Point Tools or it is contained in this toolbox as well

So, regardless of the implementation, they can be illustrative in exploring point patterns and generating containers to describe them.  I recently posted on the Minimum Area Bounding Ellipse (MABE) and I have a toolset that produces Minimum Area Bounding Circles (MABC), rectangles (MABR) and extent rectangles.  This implementation of the concave hull was proposed by Adriano Moreira and Maribel Yasmina Santos in about 2007 and there are several implementations in various languages on GitHub for those needing a trolling session.

The 'tightness' of the concave hull by changing the number of nearest neighbors to include when you are trying to decide on which points on the perimeter to keep or dump.   This 'K' factor illustrates some of the possible outcomes.  The thing to watch out for is producing degenerate points which are outside the hull, but are just to much of an outsider to be allowed into the fold.  So in pictoral form,

3 neighbours

7 neighbors

11 neighbors

And finally... with the convex hull surrounding it.

And if that isn't tight enough for you, radial sorting and concave hull generation is as tight as it is going to get

So if you need to 'contain' something, add the concave hull to your suite of tools.

I will add the implementation of concave and convex hull to my

18 Comments
MeghanBogaerts
New Contributor II

This is incredibly useful, I'm surprised that the concave hull concept hasn't been incorporated directly into Arc. Unfortunately, I don't have access to Pro. Do you have a version of this toolset available for ArcMap?

DanPatterson_Retired
MVP Emeritus

The required dependencies for ArcMap's version are not new enough.  NumPy must be greater than the distributed version of 1.9.

MeghanBogaerts
New Contributor II

Fair enough. There's a QGIS plugin that does the trick as well, I was just hoping to incorporate into a larger Arc workflow. Thanks anyway for writing up the module.

SimonJackson
Occasional Contributor III

Hi Dan.

  • I have downloaded your Point Tools and imported into my Project.
  • Running 2.4 of ArcGIS Pro
  • Trying to create a concave hull on 14k points (will later be expanding the analysis to 2M points)
  • Have tried a few different tests, but I get the following error

Traceback (most recent call last):
 File "C:\_Arcgis_Util\PointTools_pro\Scripts\hulls.py", line 331, in <module>
 p = p.view(np.float64).reshape(n, 2)
ValueError: cannot reshape array of size 5700 into shape (1900,2)
 Failed to execute (ConcaveHull).

I imagine it has an issue with the data I am feeding into it?

Using the minimum bounding geometry - convex hull works ok, but I need a tighter result, so concave looks like the right approach.

Any ideas?

DanPatterson_Retired
MVP Emeritus

I need to see what you are feeding in.  A point should only have an XY coordinate, perhaps you have a z-enabled or m-enabled point file.

Also, k=11... don't start there, start at k=3.  The higher you go, the weirder the results may be, it doesn't necessarily make the hull any tighter.

And if you have 2 million points, you would be advised to remove some duplicates or interior points if you can via selecting by distance to lines (longish story).

I would also try one of your groupings rather than hoping for the one button solves all approach since total time can be greatly shortened by doing in small batches,  It also enables you to see whether it is returning what you want since there is "No One Concave Hull", they are estimates of what could be one solution

SimonJackson
Occasional Contributor III

Thanks Dan.  Some info:

Traceback (most recent call last):
  File "C:\_Arcgis_Util\PointTools_pro\Scripts\hulls.py", line 331, in <module>
    p = p.view(np.float64).reshape(n, 2)
ValueError: cannot reshape array of size 6906 into shape (2302,2)
 Failed to execute (ConcaveHull).

"remove some duplicates or interior points if you can via selecting by distance to lines"

Good tip.  I was initially contemplating running the script in batches.

More background:

This forms part of a larger process for an accessibility workflow, that I will aim to script up and ideally make reusable for other cities.  For this initial project, looking at accessibility to schools in order to come up with some better school zones (Voronoi polygons is how they currently decide which school properties are assigned to...)

  • Property points + location of interest (my study will be schools, but could be transit stops, job centers, hospitals, etc.)
  • Origin Destination Matrix - Origins: Properties,  Destinations: Schools  (Options in here around travel mode, transit/walking/etc)
  • This then gives you the result you saw in my first screenshot - properties, nearest school ID, including travel time to the nearest school - screenshot symbolised by the school id (DestinationID) - (layer package 1.2MB)
  • Next step is to create the regions around the OD results.
    • Convex Hull works fine with the data, but the results are not suitable - overlapping
    • My initial thinking to then create the zones, is to make use of a convex hull.  Not going to work
    • Concave hull looks suitable.  
    • I would imagine there will still be some topological tests, possibly using integrate to ensure a seamless region layer at the end.  

Convex Hull - not suitable:

Let me know if you have any ideas why your tool might not be working and if you have any additional ideas I could try.  I will be sharing my workings once I hopefully get to a happy end result.

DanPatterson_Retired
MVP Emeritus

Simon... short fix since I am on Pro 2.4, beta 2 right now

In 'hulls.py'  you have to do some editing

IFFFFF you do can run this in your IDE command prompt (ie spyder's IPython consol) and line 2 doesn't fail...

import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu

then add line 2 to the script after the numpy import

Now lines 398ish edit/add to read

for i in range(0, len(groups)):
    p = groups[i]
    p = p[['SHAPE@X', 'SHAPE@Y']]
    n = len(p)
    p = stu(p)  # ---- big Stu comes to the rescue

I can remember when structured_to_unstructured was added to numpy (I am using numpy 1.16.4 right now, but I think 1.16-ish)

Let me know your version of Pro and numpy ( In [1]  numpy.version.version )

If it works out, I will update Point tools and the standalone version.

SimonJackson
Occasional Contributor III
  • I'm also running 2.4 Beta 2.
  • numpy.version.version =  '1.16.2'

My hulls.py (from your point tools toolbox) initially has this code (l.328):

for i in range(0, len(groups)):
    p = groups[i]
    n = len(p)
    p = p[['SHAPE@X', 'SHAPE@Y']]
    p = p.view(np.float64).reshape(n, 2)

Tweaking this to:

for i in range(0, len(groups)):
    p = groups[i]
    p = p[['SHAPE@X', 'SHAPE@Y']]
    n = len(p)
    p = stu(p)  # ---- big Stu comes to the rescue
  • Works for both the smaller single OD result, and for the larger test area.  Will have a look at crunching the 3M points once I have finessed some of the travel mode settings.

Side question.  I was hoping it might carry through the grouping field into the resulting polygon layer.

I can potentially get around this if I end up iterating through the OD destination results, and apply some update fields in the process.  But I imagine that would be a useful attribute to bake into the result?

Thanks for your quick responses Dan.  I think my owed beer tab for you is going to get me into debt

TKM_
by
New Contributor II

Hi Dan, could there be a way forward for Convex Hulls? #convex hull

DanPatterson_Retired
MVP Emeritus

convex hulls? numerous implementations, even mind

Point Tools convex hull is under the analysis toolset

Then esri has 

Minimum Bounding Geometry—Data Management toolbox | ArcGIS Desktop 

but for some reason, the convex hull requires an advanced license.

Keep an eye on my blog... I just put 

/blogs/dan_patterson/2019/06/26/free-tools-feature-extent-to-poly-features 

which normally requires an advanced license... I guess I should round it out with the other 'container' types.

williamwinner
New Contributor III

The problem that I've always had with this traditional implementation of a concave hull is that for most datasets that I use, it reaches the maximum k and fails.  I have never gotten it to work consistently.

I am currently using my own way to approximate a concave hull for Army Corps of Engineers survey datasets.  The data consists of thousands of points arranged in a line with lines more or less evenly spaced a few several hundred meters (sometimes up to 1000 meters) apart:

USACE sample data

So, I need to know the extents of the data.  I wrote a series of python scripts that accomplishes what I need and I have either an aggressive or a conservative approach to it.  Here are the general steps:

  1. Determine statistics on the distance between nearest points, specifically the mean and standard deviation.
  2. Then I buffer (using shapely and fiona because ESRI buffer is incredibly slow) using mean + 4 standard deviations.  This identifies lines.
  3. Then, I determine the distance (and standard deviation) between lines using the polygon centroids.
  4. From there, I re-buffer the points using mean + 1 standard deviation.  This gets me a buffer around the features that looks like this:
    USACE data buffered
  5. Then, my conservative approach is to take that buffered polygon and buffer by half the negative of the mean + 1 standard deviation.  That's not really a bounding polygon but it's close enough for my purposes.
  6. The aggressive approach is to take the second buffer and inverse buffer by the distance.
    inverse buffer
  7. That creates a lot of concave lines which I really don't want.  So, I needed to figure out a way to flatten those lines.  So, I do the following:
    1. Convert polygon to line
    2. Split the lines based on the original points and the original mean distance between points.  That more or less splits the lines at each survey line.
    3. Then, I run the Simplify Line using BEND_SIMPLIFY and the mean distance between survey lines.
    4. Then I convert that to a polygon.
    5. Originally, I thought this worked pretty well but on closer look, I noticed some points were outside so I re-buffered the that polygon by the mean distance between points.  In the end, I get this:
      Final Minimum Bounding Polygon

It may not meet everyone's needs but, for my purposes, it works great.  There are still areas like:

Imperfect?

The real question is, how should it handle something like this?  There's no one-size fits all solution so sometimes you just have to figure out a way to make it work.

DanPatterson_Retired
MVP Emeritus

Have you tried

  • group your point sequences (kd tree or by id?)
  • convert to polylines, then take the start and end points of each line sequence
  • compute the concave hull as a polyline
  • the last image shows a degenerative case of a curve at the end of a line.  You could then try doing a spatial selection of the input points with a finite distance and rerun the concave hull

As the literature often points out, there is no hull and there are many potential solutions. 

Thanks for posting.

MapDivision
New Contributor

I have tested both MBG in ArcGIS Pro and the script you have suggested in both cases the output leaves some important points which are basically way points collected to construct the boundary. In our use case, we want the script to use all the points on the periphery. We have a fields like, Unique Id, Time-Stamp, Direction (degrees) in the points attributes. Please suggest a way out.

DanPatterson_Retired
MVP Emeritus

Start a question on GeoNet and post some sample data with the options you used to generate your output.

MapDivision
New Contributor

Thank you , I will do that.

PhilLarkin1
Occasional Contributor III

Links on this post are broken. 

PhilLarkin1
Occasional Contributor III

Thanks for putting this out on Github! 

https://github.com/Dan-Patterson/Tools_for_ArcGIS_Pro

 

 

JanineF
New Contributor

Hi I would like to obtain the tightest concave hull, as in your last shown image. How can I get that? What do you mean by radial sorting? Thanks in advance

JanineF_0-1669109635449.png

 

About the Author
Retired Geomatics Instructor at Carleton University. I am a forum MVP and Moderator. Current interests focus on python-based integration in GIS. See... Py... blog, my GeoNet blog...
Labels