Wie kann ich eine flächengewichtete Summe aus einem Raster in ein Polygon in R extrahieren?

Question

Ich habe ein Raster (der Klasse "RasterLayer"), das zählt Leute, N. Ich habe Polygone administrative Grenzen (der Klasse "SpatialPolygonsDataFrame").

Ich möchte die Werte aus der Raster-Ebene in die Polygone extrahieren, und zwar in einer Weise, die durch den Anteil der Fläche des Raster-Kachels, der das Polygon überlagert, gewichtet wird. Angenommen, es gibt vier Rasterkacheln, deren Werte (4, 8, 12, 16) ein einzelnes Polygon überlappen. Wenn 25% jedes Rasters mit dem Polygon überlappt, würde ich erwarten, dass ihre Beiträge zur Polygonsumme (1, 2, 3, 4) betragen und die gewichtete Gesamtsumme, die in dieses Polygon extrahiert wird, im Wesentlichen gleich 10 ist ich möchte folgendes erreichen:

some_polygons$n_people = extract(count_raster$n_people, 
            some_polygons, 
            fun=sum, 
            weights=T, 
            na.rm=T) 

Allerdings, wenn ich diesen Befehl in R ausführen, erhalte ich folgendes:

„Warnmeldung: in .local (x, y,. ..): "Spaß" wurde auf "Mittelwert" geändert, andere Funktionen können nicht verwendet werden, wenn "Gewichte = WAHR" "

R wird zwangsweise auf fun = Mittelwert geändert, wenn die Gewichte verwendet werden. Ich muss Gewichte verwenden, um den Anteil des Rasters zu berücksichtigen, der in das Polygon fällt, aber ein Mittelwert lässt mich nicht die Menge berechnen, die ich möchte. Ich muss einen Weg finden, die Anzahl N im Raster zu addieren, gewichtet nach der Fläche im Raster, die in das Polygon fällt (eine gewichtete Summe, NICHT ein gewichteter Durchschnitt). Weiß jemand, wie ich das erreichen kann?

Ich bin glücklich zu versuchen, andere räumliche Datenklassen zu verwenden, um dies zum Funktionieren zu bringen. Ich habe versucht, das Raster als SpatialPixels, SpatialPoints und SpatialGrid darzustellen und mit der sp::over() Funktion zu arbeiten, aber konnte nicht herausfinden, wie man die gewichtete Summe berechnet, die ich benötige. Hilfe wird sehr geschätzt!

Answer 1

Zunächst bitte ein reproduzierbares Beispiel machen, würde es uns helfen, Sie

* EDIT für raster v2.3-12 Dies gilt sowohl für die Version raster v2.5-8 war aber nicht ;-) helfen!

Die weights werden wie folgt erläutert:

Falls TRUE und normalizeWeights = FALSE, die Funktion zurück, für jedes polygon, eine Matrix mit den Zellenwerten und den ungefähren Anteil von jede Zelle, die ist durch das Polygon abgedeckt (abgerundet auf 1/100). Wenn TRUE und normalizeWeights = TRUE sind, werden die Gewichtungen so normalisiert, dass sie zu eins addieren. Die Gewichte können zur Mittelwertbildung verwendet werden; siehe Beispiele. Diese Option kann nützlich sein (aber langsam), wenn die Polygone klein sind in Bezug auf die Zellen Größe des Raster * Objekt

So könnte man diese Sie gewichtete Summe zu berechnen, verwenden:

library(raster) 
library(sp) 

# Reproducible example 
set.seed(13) 
N = raster(nrows=4, ncol=4, xmn=0, xmx=4, ymn=0, ymx=4) 
N[] = runif(16, 50, 100) 

Ps1 = Polygons(list(Polygon(cbind(c(0.5,2.8,3.7,1.2,0.5), c(2,3.4,1.2,1.1,2)))), 
       ID = "a") 
SPs = SpatialPolygons(list(Ps1)) 
poly = SpatialPolygonsDataFrame(SPs, data.frame(onecol = c("one"), 
               row.names = c("a"))) 
# See plot below 
plot(N) 
plot(poly, add=T, border='red') 

# Extract with the arguments to get the appropriate weights 
# Check the version of your raster package for normalizeWeights! 
myextr = as.data.frame(extract(N, poly, weights=T, normalizeWeights=F)) 
#  value weight 
# 1 69.48172 0.16 
# 2 98.10323 0.08 
# 3 50.54667 0.61 
# 4 78.71476 0.99 
# 5 88.21990 0.17 
# 6 93.66912 0.17 
# 7 52.05317 0.87 
# 8 83.05608 0.85 
# 9 93.91854 0.43 

# compute your weighted sum 
mywsum = sum(myextr$value * myextr$weight) 
# [1] 314.9164 

* VOR EDIT

Mit raster v2.3-12, das Argument war anscheinend nicht vorhanden, daher musste ich die Aufgabe manuell erledigen, mit einer Warnung über die Auflösung des Rasters im Vergleich zur Größe des Polygons (dh eine komplett geschlossene Zelle wurde benötigt, sonst wären Anpassungen nötig gewesen erforderlich). Daher die zwei ersten Kommentare unter dieser Antwort.

Ps2 = Polygons(list(Polygon(cbind(c(0.5,2.8,3.7,1.2,0.5), c(2,3.4,1.2,0.2,2)))), 
       ID = "a") 
SPs2 = SpatialPolygons(list(Ps2)) 
poly2 = SpatialPolygonsDataFrame(SPs2, data.frame(onecol = c("one"), 
                row.names = c("a"))) 
plot(poly2, add=T, border='blue', lty=3) 

myextr2 = as.data.frame(extract(N, poly2, weights=T)) 

# compute the weight (cells fully enclosed = 1) 
myextr2$weight2 = myextr2$weight/max(myextr2$weight) 
#  value  weight weight2 
# 1 69.48172 0.027777778 0.16 
# 2 98.10323 0.013888889 0.08 
# 3 50.54667 0.105902778 0.61 
# 4 78.71476 0.171875000 0.99 
# 5 88.21990 0.029513889 0.17 
# 6 93.66912 0.052083333 0.30 
# 7 52.05317 0.173611111 1.00 
# 8 83.05608 0.173611111 1.00 
# 9 93.91854 0.090277778 0.52 
# 10 94.52795 0.003472222 0.02 
# 11 78.31402 0.107638889 0.62 
# 12 79.67737 0.048611111 0.28 
# 13 68.22573 0.001736111 0.01 

mywsum2 = sum(myextr2$value * myextr2$weight2) 
# [1] 428.2086 

Dies zeigt, dass das raster Paket schon groß war, aber immer noch verbessert :-D

[Warnung: Versuch normalizeWeights mit raster v2.3-12 zu verwenden, war es nicht abstürzt noch einen Fehler zu werfen, aber nicht das tun Job, also Vorsicht und aktualisieren Sie Ihre Version raster!]