Создайте буферы области, которые ограничены береговой линией

Я пытаюсь использовать ArcGIS 10.2 для создания буфера точек на основе заранее определенной области (например, 400 кв. Км). Кроме того, буферы некоторых точек находятся близко к береговой линии, что требует обрезки буферов на береговой линии и по-прежнему имеет ту же площадь, что и внутренние (400 кв. Км).

Кто-нибудь знает, как это можно сделать с помощью Model Builder или Arcpy?

У меня ограниченные навыки в Arcpy и R, но я был бы рад поработать над некоторыми сценариями, чтобы найти решение для этого.

Пожалуйста, смотрите изображение ниже, показывающее графическое представление того, чего я пытаюсь достичь

Площадь кругового буфера является монотонно возрастающей функцией радиуса буфера (в любом случае, в плоской системе координат). Таким образом, простая стратегия поиска может найти радиус Rтаким образом, чтобы область буфера радиуса, Rобрезанная до многоугольной области, Aбыла (с точностью до некоторого допуска) s.

Простейшим алгоритмом поиска будет просто бинарный поиск. Начните с двух радиусов, один очень маленький, а другой очень большой, чтобы область, которую вы хотите, находилась где-то между областью ограниченных буферов этих радиусов. Затем просто возьмите их среднюю точку и вычислите буферные области и выясните, находится ли желаемый радиус выше или ниже средней точки. Обновите свои пределы радиуса и повторяйте, пока не достигнете некоторого допуска от желаемой области.

Написание бинарного поиска на Python и использование ArcGIS Python API - хороший способ обучения! Я вполне уверен, что я сделал это в R, много лет назад ...

Вот немного кода R:

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))
}

Применение:

Сначала настройте простой полигональный регион Великобритании:

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)

Теперь определите точку:

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

Тогда вы просто:

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

Это список из двух компонентов, bэто буфер:

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

и у него есть правильная область:

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

и rявляется выводом из uniroot, который включает в себя значение радиуса буфера.

> b$r$root
[1] 63338.88

Таким образом, в этом случае ширина буфера была чуть меньше 64 км.

Единственное, с чем здесь нужно поиграться, это нижнее и верхнее начальные значения - я думаю, что вы можете интуитивно понять меньший радиус, так как sqrt(A/pi)верхний не так важен, так как алгоритм поиска будет увеличивать его, пока он не захватит интервал.

Алгоритм поиска может потерпеть неудачу, если начальный максимальный радиус действительно слишком велик, поскольку вы можете буферизовать всю область с огромным радиусом, и в этом случае изменение радиуса не изменит область ... Но разумные пределы должны остановить это.

Это практически невозможно из-за расположения точек. Вы можете создать буферы размером 400 км ² , но точки, расположенные ближе к береговой линии, всегда будут иметь меньшую площадь по сравнению с точками, расположенными дальше (> 400 км ² ).

Единственное, что вы можете сделать, - это выполнить анализ буферов на точках и впоследствии обрезать созданные буферы с помощью функции береговой линии.