Create area buffers that are limited by a coastline

The area of a circular buffer is a monotonically-increasing function of buffer radius (on a planar coordinate system anyway). So a simple search strategy can find a radius R such that the area of the buffer of radius R clipped to polygonal region A is (up to some tolerance) s.

The simplest search algorithm would just be a binary search. Start with two radii, one very small and one very big, such that the area you want is somewhere between the area of the clipped buffers of those radii. Then just take the midpoint of those and compute buffer areas, and figure out if the radius you want is above or below the midpoint. Update your radius limits and repeat until you get within some tolerance of your desired area.

Writing a binary search in Python and using the ArcGIS Python API sounds like a good way to learn! I'm fairly sure I've done this in R, years ago...

Here's some R code:

cropareabuff <- function(pt, region, target){
    f = function(r){
        b = rgeos::gBuffer(pt, width=r)
        return(gArea(gIntersection(b, region)) - target)
    }
    f
}

buff_with_area <- function(pt, region, target, lower, upper){
    f = cropareabuff(pt, region, target)
    r = uniroot(f, lower=lower, upper=upper, extendInt="upX")
    list(r=r, b=gIntersection(rgeos::gBuffer(pt, width=r$root), region))
}

Usage:

First set up a simple UK polygonal region:

library(raster); library(rgeos); library(rgdal)
uk = getData("GADM", country="GBR", level=0)
uk = spTransform(uk,CRS("+init=epsg:27700"))
uk = gSimplify(uk, tol=1000)

Now define a point:

p = SpatialPoints(coords=list(x=269042, y=235937), proj4string=CRS("+init=epsg:27700"))

Then you just:

b = buff_with_area(p, uk, 10000000000, 1, 10000)

This is a list with two components, b is the buffer:

plot(b$b, col=2)
plot(uk, add=TRUE)

and it has the right area:

gArea(b$b)
[1] 1e+10

and r is the output from uniroot, which includes the buffer radius value.

> b$r$root
[1] 63338.88

So in this case the buffer width was a little under 64km.

The only things to fiddle with here are the lower and upper initial values - I guess you can intuit a lower radius as sqrt(A/pi) and the upper isn't that important as the search algorithm will increase it until it captures the interval.

The search algorithm might fail if the initial max radius is really too large, since you might be buffering your entire region with a huge radius, in which case changing the radius wont change the area... But sensible limits should stop this happening.