Wednesday 30 May 2012

Crop Circles

There was a question on the GIS group on Linked In about creating circular buffers that were enclosed within a polygonal parcel yet of a given size. I reckon the way to do it is to try various values of the radius, computing the overlay, and repeating until you find the radius that gets you the area. The area is strictly increasing with radius, and as long as the polygon is big enough there will be a single unique solution between the min and max values. Simple. Here's the R code that does it:
clipCir <- function(parcel,x,y,area){
### if the parcel is smaller than the target area, we're stuffed:
  if(parea(parcel)<area){
    stop("cant be done, parcel too small")
  }
### radius is somewhere between a bit less than sqrt(a/pi)
### and the diagonal of the parcel.
  bb = bbox(parcel)
  rmin = sqrt(area/pi)*0.99
  rmax = sqrt(diff(bb[1,])^2+diff(bb[2,])^2)
### this function returns the difference between the area
### of the clipped poly and the target area:
  f = function(r){
    area-clipAreaR(parcel,x,y,r)
  }
### uniroot computes f for r between rmin and rmax until f~=0,
### which is when the area of the clipped poly = the target area
  soln = uniroot(f,c(rmin,rmax))
### return the polygon
  poly = makeClipP(parcel,x,y,soln$root)
  poly
}

clipAreaR <- function(parcel,x,y,r){
  return(parea(makeClipP(parcel,x,y,r)))
}

makeClipP <- function(parcel,x,y,r){
### create a polygon of a circle clipped to a polygon.
### the circle has 360 points round it, radius r, centre (x,y)
  theta = seq(0,2*pi,len=360)
### hoop jumping to make a SpatialPolygons object
  c = Polygon(cbind(x+r*sin(theta),y+r*cos(theta)))
  pc = Polygons(list(c),ID="A")
  spc = SpatialPolygons(list(pc))
  proj4string(spc)=proj4string(parcel)
### intersect them
  gIntersection(spc,parcel)
}

parea <- function(sp){
### compute the area of a spatial object by summing the
### area of the components
  sum(unlist(lapply(sp@polygons,function(x){x@area})))
}
The code relies on sp objects and the rgeos package. The uniroot function from the base package does the root finding.
parcels = readOGR("Data/","counties")
p = clipCir(parcels[2,],xy[1],xy[2],0.3)
plot(parcels[2,])
plot(p,add=TRUE,lwd=4)
Here's one solution:
The thin outline is the county boundary, and the thick line is a polygon centred at xy[1],xy[2] with area 0.3 (actually 0.299984). It should be easy enough to convert this to any language that can handle geometry objects - the only tricky part might be writing a uniroot function - the easiest way is to start at the min and max, and do interval bisection - compute the function at (min+max)/2 and then decide whether the zero is left or right of the midpoint. This used to be an exercise for our first-year Fortran programming students...

Here's another example with a weirder-shape parcel - here the required area was 0.005, and the returned polygon was 0.004999925 units.
Some brief timings tell me you can do about 3000 of these every minute on my 2yo home PC.

This is all good for a fixed centre, but someone mentioned the possibility that the centre isn't fixed, in which case you need some other condition for the shape you are finding - maximising some compactness parameter perhaps. Given that the parcel could be a complex, concave polygon (with holes?) this is a very hard problem. But I think I've solved the original question.

1 comment:

  1. Hello barry,
    It's really good article. but I don't want to use R-Code. so how to done this using python script for arcmap 10.

    Pl. help me on this.
    Thanks.

    ReplyDelete