Tuesday, 13 December 2016

Extracting values from raster stack and aggregating results using R



I'd like to extract values from a rasterstack, cell by cell, and create a list of the values:


ext <- extent(0,100,0,100)
r1 <- raster(nrows=100, ncols=100,ext)
r1[] <- sample(seq(from = 1, to = 6, by = 1), size = 10000, replace = TRUE)
r2 <- raster(nrows=100, ncols=100,ext)
r2[] <- sample(seq(from = 1, to = 6, by = 1), size = 10000, replace = TRUE)
r3 <- raster(nrows=100, ncols=100,ext)
r3[] <- sample(seq(from = 1, to = 6, by = 1), size = 10000, replace = TRUE)
r4 <- raster(nrows=100, ncols=100,ext)
r4[] <- sample(seq(from = 1, to = 6, by = 1), size = 10000, replace = TRUE)

s <- stack(r1,r2,r3,r4)

using extract I can create a matrix telling me the values in each layer of the stack, cell by cell:


y <- extract(s,c(1:ncell(s)))

now what I want to do is return only the unique rows of the matrix (unique(y)) with a count of how many times that row occurs.


In other words, stick a pin through every cell in a rasterstack and count how many times each 'vector' of values occurs, and also what those vectors are.



Answer



I think I've found the answer.


Convert to data frame and use count in the plyr package;



y <- extract(s,c(1:ncell(s)))
ydf <- as.data.frame(y)
head(ydf)
count(ydf, vars = c("layer.1","layer.2","layer.3","layer.4"))

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...