[GIS] Counting layer number of raster stack with value under some conditions and write it into output layer

rraster

I have generated a 8-day resolution gross primary production (GPP) fro several years using MODIS data.Each year has about 46 GPP layers. Now, I want to calculate the carbon uptake periods (in which GPPis larger than 1 g C m-2) based on these data. The result is a single raster layer with the pixel values indicate the number of the raster layers that is larger than 1. I have found a similar answer, but I still cannot solve it.Any response would be help.

The following code is one I used.

GPP2000h26v04=brick(GPP2000049,GPP2000057,GPP2000065,GPP2000073,GPP2000081,GPP2000089,GPP2000097,GPP2000105,GPP2000113,
GPP2000121,GPP2000129,GPP2000137,GPP2000145,GPP2000153,GPP2000161,GPP2000169,GPP2000177,GPP2000185,GPP2000193,GPP2000201,
GPP2000209,GPP2000217,GPP2000225,GPP2000233,GPP2000241,GPP2000249,GPP2000257,GPP2000265,GPP2000273,GPP2000281,GPP2000289,
GPP2000297,GPP2000305,GPP2000313,GPP2000321,GPP2000329,GPP2000337,GPP2000345,GPP2000353,GPP2000361)

CUP.days=calc(GPP2000h26v04,fun=function(x,na.rm) sum(x[x>1]))

Best Answer

You can use either calc or stackApply to perform this kind of operations. See the reproducible example below.

library(raster)

# Generate rasterBrick of random values
r1 <- r2 <- r3 <- raster(ncol=10, nrow=10)
r1[] <- rnorm(ncell(r1), 1, 0.2)
r2[] <- rnorm(ncell(r2), 1, 0.2)
r3[] <- rnorm(ncell(r3), 1, 0.2)
s <- brick(r1,r2,r3)

# StackApply
b1 <- stackApply(s, indices=1, fun=function(x, ...){sum(x > 1)})

# calc
b2 <- calc(s, fun=function(x){sum(x > 1)})

# Compare outputs
all.equal(b1, b2, values = TRUE)
Related Question