# 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, 2016 at 12:01
• How would you make the areas bigger when you cut then? By extending the radius of the buffer? Jan 27, 2016 at 12:13

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 = 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)
``````

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.

• 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. Jan 27, 2016 at 12:45
• The `rgeos` package and its `gBuffer` function, most likely... Jan 27, 2016 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. Jan 27, 2016 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). Jan 27, 2016 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) Jan 27, 2016 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). Jan 27, 2016 at 12:27
• Its not impossible, and its not even hard! Jan 27, 2016 at 12:39