0

I have multiple NetCDF files that contain crop and water data. These have around 30 levels of output (for 50 years). I would like to import these and sum them afterwards in R.

Currently, I am importing the levels as follows:

library(raster)

Crop1 <- Brick("Netcdfile", varname = yield, level)
Crop2 <- Brick("Netcdfile", varname = yield, level = 2)
CropN <- Brick("Netcdfile", varname = yield, level = N)

Crop1_mean <- mean(crop1[[1:50])
Crop2_mean <- mean(crop2[[1:50])
CropN_mean <- mean(cropN[[1:50])

Crop_Sum <- mosaic (Crop1_mean, Crop2_mean,  CropN_mean, fun = sum)

This is quite a lot of code to write for 30 levels when you have to do it for various projections/variables. Is there a way to write a loop for the importing of the 30 levels and sum the rasters afterwards?

ahmathelte
  • 559
  • 3
  • 15
Charles
  • 3
  • 1

1 Answers1

0

You can do this in several ways - also faster and smarter if you're going to do it for several variables. However, here is a small solution for you to get you started:

library(raster)

crop_sum <- 0

for (each_level in 1:30) {
  
  crop <- Brick("Netcdfile", varname = yield, level = each_level)
  crop_sum <- crop_sum + mean(crop[1:50])

}
harre
  • 7,081
  • 2
  • 16
  • 28
  • Maybe I am overlooking something simple, but this solution outputs crop_sum as 1 NA value not as a rasterbrick? It also only reads crop as level 30 of the levels, but I guess that is normal? – Charles Feb 28 '22 at 14:45
  • Hi Charles. I edited the answer an hour ago to fix a code typo. Please check that you've used the latest answer. If this doesn't work, I need a reproducible example and/or more info on what you need: I am not using raster bricks, so this is a more general solution that might need tweaks. – harre Feb 28 '22 at 15:42
  • Hi Charles. Sorry, I have just looked into the raster package more thoroughly, but in every situation mean(crop1[1:50]) will return a numeric, not a raster layer. Now I understand what you want in terms of the mosaic layer (instead of a normal sum), but that means that you need to take mean(crop1) instead. In other words, you'll need to explain more clearly what you need in terms of levels and years or it'll be difficult to help you. – harre Feb 28 '22 at 15:59
  • Hi Harre, I would like to extract all levels within the NetCdf file (using raster::stack or brick) preferably using a loop. Afterwards I would like to sum these raster bricks into one. Then I can use cellStats(x, sum) function in order to obtain the numeric value, but it still allows me to use other analysis on the raster brick created. Thank you for the help! – Charles Mar 01 '22 at 16:07
  • Regarding the mean(crop1[1:50]) issue I am a bit confused. Since I this function returns 1 raster layer for me that seems to be the average over 50 years. – Charles Mar 01 '22 at 16:09