I have several thousand irregularly shaped polygons in a shapefile. I want to be able to split each polygon into three areas, and to specify what the size of those areas (they sum to the previous total area). It does not matter what the sub-polygon shape is, as this is for visualization purposes.
How can I go about doing this? Is there a standard algorithm that exists that I can use?
One approach I considered was to get all the points that make up the polygon. Then, I would join two randomly together using a straight line, split the polygon, then check if the area was within a satisfactory tolerance. If it was too small, I would change the point in one direction; if it was too large, I would change to a point in the opposite direction.
Answer
This problem has many valid solutions. One of them works a little like your description, but instead of slicing the polygons at "random" locations you can do it purposefully in a way designed to minimize the amount of computation.
Here is the basic algorithm. Its input consists of any plane sweep direction, a polygon P of nonzero area, a target area a between zero and the polygon's area, and a nonnegative threshold t (in units of area). Its purpose is to split P with a line perpendicular to the sweep direction into two parts, one to the right of the line and the other to the left of the line, such that the difference between the righthand area and the target area a is no greater than t.
Let L be any oriented line perpendicular to the sweep direction. Define f(L) to be the area of P found at the right of L, minus a. In these terms the task is to find a zero of f. Because f is unlikely to be differentiable, but is continuous, use either a bisection method, the secant method, or--my favorite--Brent's method. All are simple and guaranteed to converge. Use t for the convergence tolerance for the argument.
That's it. Let's consider what goes into coding this. The root finding is routine--you can use a generic chunk of code for it--so the GIS work comes down to coding f. Doing so requires
1. Splitting the polygon by a line.
2. Computing the area of the piece(s) to the right of the line.
Both operations are implemented in almost any vector-based GIS. If not, you can replace the line by a very large rectangle representing the half-plane to the right of the line. Step 1 becomes
1'. Clip the polygon to the rectangle.
That is a really basic operation.
To get started with root finding, you need to find an interval in which the zero of f is guaranteed to lie. This is easy: project the polygon's envelope ("bounding box") into the direction of the line sweep. The projection is the interval you want.
This question has a long history. I implemented this algorithm for ArcView 3.x long ago and described it many times in the old ESRI user forums. Google
huber split polygon site:forums.esri.com
for discussions, links to code, enhancements and variations (such as splitting polygons into parts of desired sizes which are as compact as possible) and algorithms for raster data.
Here is what the continental US states look like (in an equal area projection) with the bottom third of each state shaded. Evidently the sweep direction was vertical.
No comments:
Post a Comment