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.

knitr + cactus + TwitterBootstrap + Jquery

A few notes and tips from a week of preparing some course notes.

I switched from Sweave to knitr because all the cool kids are doing it. And its better. The caching has already saved me more time than I spent switching from Sweave to knitr, which wasn't much time at all. Win

At first I was building a PDF from .Rnw sources, but then thought maybe HTML would be a better delivery platform. No need to open a PDF reader. So I figured out the .Rhtml syntax and it was done. Win.

But how to do my equations? Mathjax. Problem solved. Almost pixel-perfect equations using LaTeX syntax. Users without Javascript get the LaTeX source so its not all lost on them. Win.

Making a simple HTML file is easy enough, but I want to make a few tutorial sessions, and some other info, so I want to use a template engine. My current fave is cactus.py - you design Django templates, then write your pages to fill in sections such as 'title' or 'content' in the template, then cactus.py will build the static pages. These can then be served by a simple web server or straight from the file system. Win.

So instead of a full HTML file, I just got knitr to process a Django template chunk. However, this didn't give me the CSS for the source code syntax highlighting because there's no HEAD section in the Html. Solved by copying the knitr CSS into a file and serving it up from my template. Works fine, and serving CSS from files rather than inline in HTML is a win because now the server or client can cache it. Win.

For style and structure, I just write my templates to use styles from the Twitter Bootstrap framework. Include the CSS and JS in your template base and you have good-looking responsive pages. Sure, they look like every other Twitter Bootstrap site, but the flexibility to style them is there too. Win.

Then I figured something real cool. The R syntax highlighting wraps all functions with a span of class "functioncall". One line of jQuery can turn all those spans into hyperlinks to information from my current favourite R documentation - the R Graphical Manual. Here's some code which might get mangled by Blogger's formatting and/or my inexperience trying to stick hypertext tags into a Blogger doc:

$("span.functioncall").replaceWith(function(){return '< a href="http://rgm2.lab.nig.ac.jp/RGM2/search.php?query='+$(this).text()+'" >'+$(this).text()+'< /a >'})

Now every function links to the search for that function name. It's probably hard to link directly to the help for that function exactly, since the system would have to know what package each function was in to get to exactly the right help on the RGM, but this is a big help. Win.

So that was my weekend. Also, Chelsea FC. Win.