Obszar okrągłego bufora jest monotonicznie rosnącą funkcją promienia bufora (w każdym razie na płaskim układzie współrzędnych). Tak więc prosta strategia wyszukiwania może znaleźć taki promień R
, że obszar bufora promienia R
przyciętego do regionu wielokąta A
wynosi (do pewnej tolerancji) s
.
Najprostszym algorytmem wyszukiwania byłoby po prostu wyszukiwanie binarne. Zacznij od dwóch promieni, jednego bardzo małego i jednego bardzo dużego, tak aby żądany obszar znajdował się między obszarem przyciętych buforów tych promieni. Następnie po prostu weź środkowy punkt tych i oblicz obszary buforowe i dowiedz się, czy żądany promień jest powyżej czy poniżej środkowego punktu. Zaktualizuj granice promienia i powtarzaj, aż osiągniesz tolerancję żądanego obszaru.
Pisanie binarnego wyszukiwania w Pythonie i używanie ArcGIS Python API wydaje się dobrym sposobem na naukę! Jestem całkiem pewien, że zrobiłem to w R, lata temu ...
Oto kod 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))
}
Stosowanie:
Najpierw skonfiguruj prosty wieloboczny region Wielkiej Brytanii:
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)
Teraz zdefiniuj punkt:
p = SpatialPoints(coords=list(x=269042, y=235937), proj4string=CRS("+init=epsg:27700"))
Więc po prostu:
b = buff_with_area(p, uk, 10000000000, 1, 10000)
To jest lista z dwoma składnikami, b
to bufor:
plot(b$b, col=2)
plot(uk, add=TRUE)
i ma odpowiedni obszar:
gArea(b$b)
[1] 1e+10
i r
jest wyjściem uniroot
, które zawiera wartość promienia bufora.
> b$r$root
[1] 63338.88
W tym przypadku szerokość bufora wynosiła nieco poniżej 64 km.
Jedynymi rzeczami, które można tutaj manipulować, są dolne i górne wartości początkowe - myślę, że można intuicyjnie ustalić niższy promień, ponieważ sqrt(A/pi)
górna nie jest tak ważna, ponieważ algorytm wyszukiwania zwiększy ją, dopóki nie zarejestruje interwału.
Algorytm wyszukiwania może zawieść, jeśli początkowy maksymalny promień jest naprawdę zbyt duży, ponieważ możesz buforować cały region z dużym promieniem, w którym to przypadku zmiana promienia nie zmieni obszaru ... Ale rozsądne granice powinny temu zapobiec.