[R-sig-Geo] Creating random (but believable) geographies

[Jim Lemon]

On 12/28/2009 09:33 PM, Barry Rowlingson wrote:
> I've recently been working on a web site for tropical disease mapping
> in Africa. As a demo, I'm not using real data, but I have been using
> real countries. I don't really want to use real countries just in case
> there's any misunderstanding that this is simulated case numbers.
>
>   So I've thought about constructing random polygons to represent a
> fake country and its administrative subdivisions. Any ideas how to do
> this? Here's some thoughts so far:
>
>   * Generate some random points, create voronoi polygons, then wiggle
> the edges to make them look more natural. Problem may be stopping the
> wiggling from causing edges to overlap, but maybe that just needs a
> bunch of point-in-poly tests... Hmmm...
>
>   * Overlay two existing unassociated polygon geometries, such as UK
> counties and Ethiopian districts, compute intersections and rebuild
> the polygon topology. Then you should end up with polygonal regions
> with outlines composed of both sets of polygons. You might end up with
> lots of tiny polygons where two lines closely overlap along their
> length, but maybe a post-processing step can fix this.
>
>   * Raster-based methods: similar to the voronoi method, start with a
> number of seed points on a grid and grow by some random method - then
> do raster-vector conversion.
>
>   Obviously there's lots of possible options such as creating islands
> and holes, but I thought I'd start with the basics.
>
>   If anyone has any existing fake geographies in digital form then it
> might be an idea to create a repository of them...
>
>    
Hi Barry,
The first thing I thought of was to specify nodes (intersections of 
borders = vertices of polygons) and then wiggle the lines between the 
nodes using something like this:

wigglyLine<-function(xstart,ystart,xend,yend,nseg=20,mult=1) {
  xvec<-seq(xstart,xend,length.out=nseg+1)
  yvec<-seq(ystart,yend,length.out=nseg+1)
  xvec[2:nseg]<-xvec[2:nseg]+(runif(nseg-1)-0.5)*mult
  yvec[2:nseg]<-yvec[2:nseg]+(runif(nseg-1)-0.5)*mult
  return(list(x=xvec,y=yvec))
}

Obviously there are a lot of things that you could fiddle with here, 
including the easiest way to specify the vertices and then turn the 
wiggled edges into polygons. It does solve the problem of having the 
polygons fit together, though.

Jim