Create area buffers that are limited by a coastline

Question

I'm trying to use ArcGIS 10.2 to create a buffer of points based on a pre-defined area (e.g 400 sq km). Further to that, the buffers of some of the points are close to the coastline that requires the buffers to be clipped at the coastline and still have the same area as the ones that are inland (400 sq km).

Does anyone know how this could be done with either Model Builder or Arcpy?

I have limited skills with Arcpy and R but would be happy to work on some scripting to get a solution for this.

Please see image below showing a graphical representation of what I am trying to achieve

Answer 1

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.

Answer 2

It's almost impossible, due to the position of the points. You can create buffers of 400km², but points closer to the coastline will always have a smaller area compared to the ones further away (>400km²).

The only thing you can do is do perform a buffer analysis on the points and clip the created buffers with the coastline feature afterwards.