I have a raster stack created from daily climate data. Can be found here:
#!/bin/bash
wget -nc -c -nd http://northwestknowledge.net/metdata/data/tmmx_1982.nc
The goal is to get the 95th percentile of temperature values per month from these daily records. Whenever I use calc from the raster package it just returns one layer instead of 12 (e.g., 12 months) What am I missing?!
library(raster)
library(tidyverse)
library(lubridate)
file = "../data/raw/climate/tmmx_1982.nc "
rstr <- raster(file)
> rstr
class : RasterBrick
dimensions : 585, 1386, 810810, 366 (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667 (x, y)
extent : -124.793, -67.043, 25.04186, 49.41686 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
data source : in memory
names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ...
min values : 1.3268673, 0.7221603, 1.8519223, 1.6214808, 0.8629752, 1.1126643, 1.8769895, 0.9587604, 1.7360761, 2.1099827, 2.1147265, 1.8696048, 1.7619936, 2.0253942, 2.6840794, ...
max values : 73.20462, 60.35675, 64.68890, 53.11994, 60.15675, 55.91125, 77.29095, 64.39179, 48.26004, 64.70559, 79.85970, 62.31242, 53.89768, 52.15949, 80.23198, ...
date_seq <- date_seq[1:nlayers(rstr)]
month_seq <- month(date_seq)
mean_tmp <- stackApply(rstr, month_seq, fun = mean)
> mean_tmp
class : RasterBrick
dimensions : 585, 1386, 810810, 12 (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667 (x, y)
extent : -124.793, -67.043, 25.04186, 49.41686 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
data source : /tmp/RtmpYf4pQe/raster/r_tmp_2017-09-25_182536_48012_88372.grd
names : index_1, index_2, index_3, index_4, index_5, index_6, index_7, index_8, index_9, index_10, index_11, index_12
min values : 4.586111, 5.656802, 6.444234, 6.875973, 6.281896, 4.495534, 5.081545, 4.396824, 4.316368, 6.413400, 4.233641, 3.119827
max values : 49.12178, 47.61632, 44.70796, 47.57829, 46.97714, 51.61986, 37.77228, 51.30043, 42.51572, 36.86453, 37.96615, 52.15552
mean_90thtmp <- calc(mean_tmp, forceapply = TRUE,
fun = function(x) {quantile(x, probs = 0.90, na.rm = TRUE) })
> mean_90thtmp
class : RasterLayer
dimensions : 585, 1386, 810810 (nrow, ncol, ncell)
resolution : 0.04166667, 0.04166667 (x, y)
extent : -124.793, -67.043, 25.04186, 49.41686 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
data source : in memory
names : layer
values : 8.84197, 50.52144 (min, max)
Suggestions are very much appreciated!
Thanks!
One option is using a for loop:
x <- stack() # create an empty stack
for (i in 1:nlayers(mean_tmp)){
mean_90thtmp <- calc(mean_tmp[[i]], forceapply = TRUE,
fun = function(x) {quantile(x, probs = 0.90, na.rm = TRUE) })
x <- stack(x , mean_90thtmp )
}