[R-sig-Geo] How to buffer a polygon by area

[Barry Rowlingson]

On Wed, May 28, 2014 at 11:37 AM, Julie Lee-Yaw <julleeyaw at yahoo.ca> wrote:
> Hello,
> Can anyone tell me how to buffer a polygon by area?
> I have a simple (no-holes and only one feature) polygon shapefile in an equal area projection and I need to buffer it to extend the boundaries so that the resulting polygon is 10% bigger by area. I know how to buffer by a fixed distance using the rgeos package but I will have to repeat this process across many different shapefiles that are of different sizes and a fixed buffer does not make sense in my case.
> Thanks in advance for the help!

 Since the area of a buffered feature must be increasing with
increasing buffer width I think you could use a binary search
algorithm on the buffer width. This might require doing five or six
buffers on each feature until you get to your required accuracy.

Here's something:

abuff <- function(f, percent, precision, maxsteps = 20){
    require(rgeos)
    A0 = gArea(f)

    targetA = A0 + A0 * (percent/100)
    e = bbox(f)
    w0 = 0
    w2 = diff(e[1,])/10 # initial w

    ## initial expansion until we pass the solution:

    repeat{
        b1 = gBuffer(f, width=w2)
        if(gArea(b1) < targetA){
            w2 = w2 * 2
        }else{
            break
        }
    }
    w1 = (w0 + w2)/2

    ## solution is between w0 and w2

    steps = 1

    repeat{
        b = gBuffer(f, width=w1)
        A = gArea(b)
        if(abs(100*(A-targetA)/targetA) < precision){
            return(list(steps = steps, poly = b, target=targetA))
        }
        if(A < targetA){
            w0 = w1
            w1 = (w1+w2)/2
        }
        if(A > targetA){
            w2 = w1
            w1 = (w2+w0)/2
        }
        steps = steps + 1
        if(steps > maxsteps){
            return(NA)
        }
    }
}

Sample usage is, given 'london' as a single row spatial polygons data frame:

 > lbuff = abuff(london, 10, .5)

 - for a 10% area increase to a 0.5% precision

 > lbuff$steps
[1] 6

 - it took 6 steps, so had to compute 6 buffers.

 > plot(lbuff$poly)
 > plot(london,add=TRUE)

 - plots the buffer and the source.

 > lbuff$target
[1] 1766883778
 > gArea(lbuff$poly)
[1] 1768395765

The target element is the area we were trying to find, and the next
line computes the area of the result. If you want more accuracy,
increase the precision at the expense of needing to do more buffers.

I've put this function in a gist if anyone wants to improve it:

https://gist.github.com/barryrowlingson/d1d21c537777cfda8178

Barry