0

I have a RasterBrick that contains mean values one layer for each month in a 72-year interval. I want to get the mean value for each year - i.e. return a 72-layer RasterBrick.

The following code has worked on other similar rasters, yielding the expected results (found here):

data <- raster::brick(".../air.mon.mean.nc", varname = "air")

index <- format(as.Date(raster::getZ(data), format = "%Y-%m-%d %H:%M:%S"), format = "%Y")
index <- as.numeric(index)

yearly <- raster::stackApply(data, index, fun = mean) 

> yearly
class      : RasterBrick 
dimensions : 360, 720, 259200, 72  (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5  (x, y)
extent     : 0, 360, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : C:/Users/villar/AppData/Local/Temp/RtmpAbUQQT/raster/r_tmp_2019-08-05_102157_18368_64365.grd 
names      : index_1948, index_1949, index_1950, index_1951, index_1952, index_1953, index_1954, index_1955, index_1956, index_1957, index_1958, index_1959, index_1960, index_1961, index_1962, ... 
min values :         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA, ... 
max values :         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA,         NA, ...

However, when running on this data, it returns only NA values.

fun = function(x, na.rm) {sum(x)/12} does not work, nor does adding na.rm = TRUE.

Any help would be greatly appreciated!

The data is downloaded from here (air.mon.mean.nc).

VLarsen
  • 67
  • 1
  • 7
  • On my system, it looks like you are making a small mistake in formatting the times: `raster::getZ(data)` returns dates in the format `"%Y-%m-%d"`, not `"%Y-%m-%d %H:%M:%S"`. This means that `index` is all `NA` and thus `stackApply` does not work properly and returns only one layer with all `NA`. In any case, this does not happen to you, and I can reproduce your results if I use `"%Y-%m-%d"`. – AF7 Aug 05 '19 at 19:32
  • For some reason, I can get it working if I do `fun = function(x, na.rm=TRUE) {sum(x, na.rm = na.rm)/12}` and then `yearly <- raster::stackApply(data, index, fun = fun)`... but then it's all 0s. – AF7 Aug 05 '19 at 19:45

2 Answers2

1

I have no idea why raster behaves like this. I can reproduce it. However, if you can use CDO, you could use the operator yearmean to get the output you want: cdo yearmean input.grb output.grb. This will also be faster than any R implementation using a single core.

If you prefer to stay inside R, I suggest you take a look at the new-but-awesome stars package.

You could do something like this:

library(stars)

s = read_stars(s)

yrs = st_get_dimension_values(s, 'time') %>% format('%Y') %>% as.numeric
mymean = function(v, indices, fun = mean, na.rm = FALSE) {
    sapply(unique(indices), function(i) fun(v[indices == i], na.rm = na.rm))
}

yearly = st_apply(s, 1:2, mymean, indices = yrs, na.rm = TRUE)

Optionally you could also do this multicore (see ?st_apply) but I doubt it would be faster than barebones CDO anyway.

AF7
  • 3,160
  • 28
  • 63
0

I had the same problem but this It worked for me calc(data,fun=function(x) { by(x, index, sum)}). Could you solve it with the stackApply function?

tmsppc
  • 107
  • 3
  • 10