Estimating the boundary of arbitrarily distributed data

Question

I have two dimensional discrete spatial data. I would like to make an approximation of the spatial boundaries of this data so that I can produce a plot with another dataset on top of it.

Ideally, this would be an ordered set of (x,y) points that matplotlib can plot with the plt.Polygon() patch.

My initial attempt is very inelegant: I place a fine grid over the data, and where data is found in a cell, a square matplotlib patch is created of that cell. The resolution of the boundary thus depends on the sampling frequency of the grid. Here is an example, where the grey region are the cells containing data, black where no data exists.

OK, problem solved - why am I still here? Well.... I'd like a more "elegant" solution, or at least one that is faster (ie. I don't want to get on with "real" work, I'd like to have some fun with this!). The best way I can think of is a ray-tracing approach - eg:

1. from xmin to xmax, at y=ymin, check if data boundary crossed in intervals dx
2. y=ymin+dy, do 1
3. do 1-2, but now sample in y

An alternative is defining a centre, and sampling in r-theta space - ie radial spokes in dtheta increments.

Both would produce a set of (x,y) points, but then how do I order/link neighbouring points them to create the boundary?

A nearest neighbour approach is not appropriate as, for example (to borrow from Geography), an isthmus (think of Panama connecting N&S America) could then close off and isolate regions. This also might not deal very well with the holes seen in the data, which I would like to represent as a different plt.Polygon.

The solution perhaps comes from solving an area maximisation problem. For a set of points defining the data limits, what is the maximum contiguous area contained within those points To form the enclosed area, what are the neighbouring points for the nth point? How will the holes be treated in this scheme - is this erring into topology now?

Apologies, much of this is me thinking out loud. I'd be grateful for some hints, suggestions or solutions. I suspect this is an oft-studied problem with many solution techniques, but I'm looking for something simple to code and quick to run... I guess everyone is, really!

Cheers,

David

~~~~~~~~~~~~~~~~~~~~~~~~~

OK, here's attempt #2 using Mark's idea of convex hulls:

For this I used qconvex from the qhull package, getting it to return the extreme vertices. For those interested:

cat [data] | qconvex Fx > out

The sampling of the perimeter seems quite low, and although I haven't played much with the settings, I'm not convinced I can improve the fidelity.

Answer 1

I think what you are looking for is the Convex Hull of the data That will give a set of points that if connected will mean that all your points are on or inside the connected points

Answer 2

I may have mixed something, but what's the motivation for simply not determining the maximum and minimum x and y level? Unless you have an enormous amount of data you could simply iterate through your points determining minimum and maximum levels fairly quickly.

This isn't the most efficient example, but if your data set is small this won't be particularly slow:

import random
data = [(random.randint(-100, 100), random.randint(-100, 100)) for i in range(1000)]

x_min = min([point[0] for point in data])
x_max = max([point[0] for point in data])

y_min = min([point[1] for point in data])
y_max = max([point[1] for point in data])