The other day I asked a question on the R-sig-geo mailing list. How can I generate a set of polygons that look a bit like a real set of districts or other subdivisions? Here's what I came up with:
- Take a set of points
- Generate the voronoi tiles of the points
- Turn the straight-line edges into something wiggly
- Clip or cut the resulting polygons
So I wrote some R code to do all this. 1 and 2 are pretty easy using the tripack package to generate the voronoi tiles. Making the edges wiggly is a bit trickier. I implemented a fractal algorithm. Split the edge into two pieces by making a new point somewhere near the centre of the edge. Then repeat for the two pieces to make four pieces. Do this three times. This should make a fairly straight wiggly line. For clipping and cutting the polygons I used the gpclib package (free for non-commercial use only - read the license).
Clipping and cutting is necessary because voronoi tiles can extend way beyond the initial points if the tile edges are nearly parallel. I implemented four clipping modes:
- No clipping (so you get the full polygons)
- Include polygons that overlap a given boundary polygon
- Clip polygons to given boundary polygon
- Include only polygons completely inside a boundary polygon
By default the boundary polygon is the convex hull of the input points, but you can specify another polygon for clipping.
Here's my sample data - some points and a blue polygon for clipping:
Here's a fake geography created using these points and clipmode zero - the blue polygon is ignored here:
Note that the polygons bottom left head off some distance south-west. The lack of polygons in the north and west is due to these areas being infinite. They don't form closed tiles. This means the number of polygons you get returned wont be the same as the number of input points.
The wiggling algorithm here may cause lines to cross and overlap. There's no test for that at the moment.
Clipmode 1 only includes polygons that are partly or wholly inside the clip polygon:
This has eliminated all the extending polygons in the south-west. Mode 2 clipping results in the polygons being clipped to the boundary polygon. This is useful if you are creating a set of polygon districts within an existing area such as a country boundary:
Note you can't see the red edges under the blue polygon boundary here - they are there though.
Mode 3 only returns polygons wholly inside the clip polygon:
This only leaves 8 polygons - it looks like one polygon on the inner corner has just been clipped out.
So there it is. I've not polished up the code yet, but it's pretty short and sweet. Interested? Let me know.
Addendum: I promise I'll bundle the code up shortly, but today I've been making tweaks to produce fake river networks:
It's pretty much the same basic procedure - generate a voronoi tiling and wiggle the edges. Then compute a minimum spanning tree (any tree will probably do) with distance as the weight. The resulting map looks a bit like a river network (although I wouldn't want to draw the contours or DEM it came from).
I tried various weighting schemes for the MST. Using distance means it uses up the small edges before the large edges, which tends to leave those big gaps that could be mountain ridges.
Other ideas welcome!
