# Create area buffers that are limited by a coastline

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 • Would you be able to include a picture of what you are trying to describe in words? – PolyGeo Jan 27 '16 at 12:01
• How would you make the areas bigger when you cut then? By extending the radius of the buffer? – Peter Horsbøll Møller Jan 27 '16 at 12:13

## 2 Answers

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)
 1e+10
``````

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

``````> b\$r\$root
 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.

• How did you do this in R? I forgot to mention that I have some experience in R so I wouldn't mind a solution using R as well. – Funkeh-Monkeh Jan 27 '16 at 12:45
• The `rgeos` package and its `gBuffer` function, most likely... – Spacedman Jan 27 '16 at 12:49
• Actually I tell a lie, I implemented something like it in Python as a QGIS plugin - it buffered polygons until the buffered poly was 2x (or Nx) the area of the original polygon. Same search algorithm though. – Spacedman Jan 27 '16 at 13:01
• +1. The advantages of the approach shown in the `R` code are (a) it separates the GIS calculations from the search logic and (b) it capitalizes on search algorithms (in `uniroot`) that have been optimized and tested--you don't have to write one yourself (and it will likely not be the most efficient). – whuber Jan 27 '16 at 15:56
• I suspect scipy implements similar root-finding algorithms in its optimize module: docs.scipy.org/doc/scipy/reference/optimize.html (yes, ?uniroot cites Brent, scipy has Brent-ish functions) – Spacedman Jan 27 '16 at 16:52

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

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.

• It may not be impossible, but it could be an NP Complete problem that may confound solution. Getting the area perfect is the challenge (may take scores of iterations to get close). – Vince Jan 27 '16 at 12:27
• Its not impossible, and its not even hard! – Spacedman Jan 27 '16 at 12:39