zonale Statistiken R

Question

Ich verwende das raster Paket in R und führen ein extract auf ein Raster mit shapefile der mehr polygons (200) enthält. Jedes Polygon hat eine Spalte, in der class des Polygons angegeben ist. Wenn ich die extract führe, erhalte ich eine dataframe, wobei jeder Wert des Rasters einem Polygon "zugewiesen" wird: In welchem ​​Polygon befindet sich welches Pixel. Die Informationen, die mein Polygon enthält (d. H. Welche Klassen die sind), sind jedoch weg. Mein Datenrahmen nach den extract sieht ungefähr so ​​aus:

test <- extract(myRaster, myTrainingPolygon, df = TRUE) 

ID  band1  band2 band3 
1  0.101  0.827 ... 
... ...  ...  ... 
200 0.876  0.821 ... 

Was ich brauche, ist so etwas wie

ID  band1  band2 band3 class 
1  0.101  0.827 ...  class1 
... ...  ...  ...  ... 
200 0.876  ...  ...  class3 

Wie erhalte ich die Informationen dort - oder in erster Linie - sie nicht verlieren, wenn Ich führe den Auszug durch ?!

Answer 1

Immer umfassen und Beispieldaten

library(raster) 
r <- raster(ncol=36, nrow=18) 
r[] <- 1:ncell(r) 
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)) 
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)) 
p <- spPolygons(cds1, cds2) 
p <- SpatialPolygonsDataFrame(p, data.frame(class=c('A', 'B'), stringsAsFactors=FALSE)) 

Jetzt verwenden wir tun können:

v <- extract(r, p) 
str(v) 

#List of 2 
# $ : num [1:38] 326 327 328 329 330 331 332 333 334 335 ... 
# $ : num [1:25] 172 173 208 209 244 245 279 280 281 282 ... 

erstellen data.frame, die das Polygon sequenzielle ID und die extrahierten Werte hat

d <- data.frame(id=rep(1:length(v), sapply(v, length)), value=unlist(v)) 

(Dieser spezielle Ansatz mit unlist funktioniert nur für eine einzelne Ebene). Für mehrere Schichten

tun

#d <- data.frame(id=rep(1:length(v), sapply(v, length)), do.call(rbind, v)) 

erstellen data.frame mit der sequentiellen ID und mit dem anderen Polygon-Attribute

pd <- cbind(id=1:length(p), data.frame(p)) 

die beide

m <- merge(pd, d) 

Answer 2

?extract zeigt an, dass das Argument sp steuert, ob die extrahierten Werte zum Datenrahmen des Geo-Objekts hinzugefügt werden sollen oder nicht. Die Angabe sp=TRUE sollte den Trick machen. Z.B .:

library(maptools) 
library(raster) 
data(wrld_simpl) 
r <- raster(extent(-180, 180, -90, 90), res=10) 
r[] <- runif(ncell(r)) 
wrld_simpl_new <- extract(r, wrld_simpl, fun=mean, sp=TRUE) 

head(wrld_simpl_new) 

##  FIPS ISO2 ISO3 UN    NAME AREA POP2005 REGION SUBREGION  LON  LAT  layer 
## ATG AC AG ATG 28 Antigua and Barbuda  44 83039  19  29 -61.783 17.078 0.9142067 
## DZA AG DZ DZA 12    Algeria 238174 32854159  2  15 2.632 28.163 0.1774097 
## AZE AJ AZ AZE 31   Azerbaijan 8260 8352021 142  145 47.395 40.430 0.3098710 
## ALB AL AL ALB 8    Albania 2740 3153731 150  39 20.068 41.143 0.3746480 
## ARM AM AM ARM 51    Armenia 2820 3017661 142  145 44.563 40.534 0.3494729 
## AGO AO AO AGO 24    Angola 124670 16095214  2  17 17.544 -12.296 0.2873931 

Die zusätzliche Spalte hat den Namen "Schicht", denn das ist der (Standard) Name der RasterLayer war.