ArcGIS 10 – Create Weighted Thiessen Polygons with VBA and VB.NET

Question

I have a point shapefile and I create Thiessen (Voronoi) polygons programmatically using this scripting syntax:

CreateThiessenPolygons_analysis (in_features, out_feature_class, fields_to_copy) 

However, each point is related with an area (i.e. the preferred size of each polygon) and I wish the Thiessen polygons to be weighted based on this field.

Is that possible and how?

Is there any relevant code in VBA?

Answer 1

There are many ways to weight distances for constructing Thiessen polygons. The basic idea in constructing them is based on comparing the distance between an arbitrary point x and two fixed points p and q; you need to decide whether x is "closer" to p than to q or not. To this end--at least conceptually--we consider the distances dp = d(x, p) and dq= d(x, q). The weighting usually occurs in two ways: the points can be given positive numeric weights wp and wq and the distances themselves can be transformed.

To make sense, the transformation (which I will write as f) should increase as the distances increase; that is, f(d') > f(d) whenever d' > d >= 0. Examples of such transformations are f(d) = d+1, f(d) = d^2 (Reilly's Law of Retail Gravitation), f(d) = 1 - 1/d (assuming all distances are less than 1), f(d) = log(d), f(d) = exp(d)-1.

We would then say x is "closer" to p than to q exactly when

f(d(x, p)) / wp < f(d(x, q)) / wq.

Notice the division by the weights, rather than multiplication: this means large weights will tend to "pull in" points at larger distances. You will see this in the running example below.

Here's the beautiful thing, and the whole point of this somewhat abstract exposition: although the resulting Thiessen regions can have complex, extremely difficult to calculate boundaries, they are relatively easy to compute using a grid-based representation. Here's the recipe:

1. For each input point p, compute its Euclidean distance grid [d(p)].

2. Use Map Algebra to apply f and the weights, thereby re-expressing each distance grid as

   [fp] = f([d(p)]) / wp.

   Here is an example using f(d) = 100 + d^(3/2); the scale is 400 by 600.

   

   As f(d) increases the value gets darker. Evidently the distance in this example is with respect to the central red point; the other four points get their separate distance calculations (not shown). The areas of the dots are proportional to their weights, which are 2, 10, 3, 4, and 5.

3. Compute the local minimum of all these grids [fp]. Call this [f]. Here is an example.

   

4. By comparing [f] to each [fp], to each grid cell assign the identifier of the first p for which [f] >= [fp]. (This can be done in one step with a lowest position operation, for instance.)

   

   (I doubt there exists an algorithm anywhere that will compute a vector-format solution for this weighting function f.)

Obviously if you have more than a handful of points p you will script this, and if their number runs into thousands you will probably abandon the attempt as being computationally impracticable (although there are ways to expedite the calculation by tiling it).

Another example, showing Thiessen polygons on an ellipsoid, appears at https://gis.stackexchange.com/a/17377/.