2

I am calculating the Maximum GPP of a year based on a 8-day GPP raster layers. I want to know how to calculate the maximum of the stacked raster at the position of each pixel and to output a new raster. I know Cell statistics in Arcmap can solve this problem,but I have to calculate a lot of data. max() and maxValue() in "raster" package are not worked.

My code is :

library(raster)
GPP2000049=raster("GPP.2000049.h23v03.tif") GPP2000057=raster("GPP.2000057.h23v03.tif")
GPP2000065=raster("GPP.2000065.h23v03.tif") GPP2000073=raster("GPP.2000073.h23v03.tif")
GPP2000081=raster("GPP.2000081.h23v03.tif") GPP2000089=raster("GPP.2000089.h23v03.tif")
GPP2000097=raster("GPP.2000097.h23v03.tif") GPP2000105=raster("GPP.2000105.h23v03.tif")
GPP2000113=raster("GPP.2000113.h23v03.tif") GPP2000121=raster("GPP.2000121.h23v03.tif")
GPP2000129=raster("GPP.2000129.h23v03.tif") GPP2000137=raster("GPP.2000137.h23v03.tif")
GPP2000145=raster("GPP.2000145.h23v03.tif") GPP2000153=raster("GPP.2000153.h23v03.tif")
GPP2000161=raster("GPP.2000161.h23v03.tif") GPP2000169=raster("GPP.2000169.h23v03.tif")
GPP2000177=raster("GPP.2000177.h23v03.tif") GPP2000185=raster("GPP.2000185.h23v03.tif")
GPP2000193=raster("GPP.2000193.h23v03.tif") GPP2000201=raster("GPP.2000201.h23v03.tif")
GPP2000209=raster("GPP.2000209.h23v03.tif") GPP2000217=raster("GPP.2000217.h23v03.tif")
GPP2000225=raster("GPP.2000225.h23v03.tif") GPP2000233=raster("GPP.2000233.h23v03.tif")
GPP2000241=raster("GPP.2000241.h23v03.tif") GPP2000249=raster("GPP.2000249.h23v03.tif")
GPP2000257=raster("GPP.2000257.h23v03.tif") GPP2000265=raster("GPP.2000265.h23v03.tif")
GPP2000273=raster("GPP.2000273.h23v03.tif") GPP2000281=raster("GPP.2000281.h23v03.tif")
GPP2000289=raster("GPP.2000289.h23v03.tif") GPP2000297=raster("GPP.2000297.h23v03.tif")
GPP2000305=raster("GPP.2000305.h23v03.tif") GPP2000313=raster("GPP.2000313.h23v03.tif")
GPP2000321=raster("GPP.2000321.h23v03.tif") GPP2000329=raster("GPP.2000329.h23v03.tif")
GPP2000337=raster("GPP.2000337.h23v03.tif") GPP2000345=raster("GPP.2000345.h23v03.tif")
GPP2000353=raster("GPP.2000353.h23v03.tif") GPP2000361=raster("GPP.2000361.h23v03.tif")

MaxGPP2000h23v03=stack(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)

MaxGPP2000h23v03=max(MaxGPP2000h23v03)
writeRaster(MaxGPP2000h23v03,filename = "MaxGPP2000h23v03.tif",overwrite=TRUE)
4

You can either use stackApply or calc to reduce a rasterStack or RasterBrick to the maximum value of each pixel.

library(raster)

# Create a test RasterBrick
r <- raster(ncol=10, nrow=10)
r[]=1:ncell(r)
b <- brick(r,r,r,r,r,r)

stackApply(b, indices = rep(1, nlayers(b)), fun = max)

calc(b, function(x){max(x)})
  • clusterR can help speed things up for larger files too. Also, for each of calc, stackApply and clusterR, you can make things still more efficient by specifying the output datatype, e.g. if you know your outputs will be whole numbers, its a good idea to specify datatype='INT2S'. The raster package creates temp gri/grd files during these operations even when you don't tell it to write to file, so its a good idea to minimise their size where appropriate. – obrl_soil Jan 12 '17 at 11:25
  • Thank you so much! It works now. I will have a look at clusterR. – dream Jan 13 '17 at 1:42
2

If your stack have NA values (or not) use:

library(raster)

maxStack <- calc(MaxGPP2000h23v03, function(x) max(x, na.rm = TRUE))

or

maxStack <- max(MaxGPP2000h23v03, na.rm = TRUE)

or

maxStack <- stackApply(MaxGPP2000h23v03, indices=rep(1,nlayers(MaxGPP2000h23v03)), fun = max, na.rm=TRUE)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.