diff --git a/DESCRIPTION b/DESCRIPTION index f4ad9b9fab9e2f60d72889f6f144952f571e3e05..dd958d8d11d54ae6df68a652db4131edfdd96860 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -41,8 +41,9 @@ Imports: Suggests: zeallot, testthat, - R.rsp -VignetteBuilder: R.rsp + knitr, + rmarkdown +VignetteBuilder: knitr License: Apache License 2.0 | file LICENSE Encoding: UTF-8 LazyData: true diff --git a/vignettes/MultiModelSkill_vignette.md b/vignettes/MultiModelSkill_vignette.Rmd similarity index 99% rename from vignettes/MultiModelSkill_vignette.md rename to vignettes/MultiModelSkill_vignette.Rmd index 6fc676109908add0004f31bdcfe2ab141b33cd4a..414955a20ff1fe8671460acd87fd629f139e5c63 100644 --- a/vignettes/MultiModelSkill_vignette.md +++ b/vignettes/MultiModelSkill_vignette.Rmd @@ -1,9 +1,9 @@ --- author: "Nuria Perez" date: "`r Sys.Date()`" -output: html_document +output: rmarkdown::html_vignette vignette: > - %\VignetteEngine{R.rsp::md} + %\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{Multi-model Skill Assessment} %\usepackage[utf8]{inputenc} --- @@ -21,7 +21,6 @@ A version of s2dverification under development should be loaded by running: devtools::install_git('https://earth.bsc.es/gitlab/es/s2dverification', branch = 'mostlikely') library(s2dverification) -library(zeallot) ``` @@ -70,8 +69,11 @@ obs <- "erainterim" Using the `CST_Load` function, the data available in our data store can be loaded. The following lines, shows how this function can be used. However, the data is loaded from a previous saved `.RData` file: Ask nuria.perez at bsc.es to achieve the data to run the recipe. -```r +```r +require(zeallot) + glosea5 <- '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc' + c(exp, obs) %<-% CST_Load(var = clim_var, exp = list(list(name = 'glosea5', path = glosea5), list(name = 'ecmwf/system4_m1'), diff --git a/vignettes/MultivarRMSE_vignette.md b/vignettes/MultivarRMSE_vignette.Rmd similarity index 98% rename from vignettes/MultivarRMSE_vignette.md rename to vignettes/MultivarRMSE_vignette.Rmd index 78af60807aa75815237fed804619e25f172de392..3712075dc8fc07c4f684bd6ea75c3367618946d2 100644 --- a/vignettes/MultivarRMSE_vignette.md +++ b/vignettes/MultivarRMSE_vignette.Rmd @@ -1,9 +1,9 @@ --- author: "Deborah Verfaillie" date: "`r Sys.Date()`" -output: html_document +output: rmarkdown::html_vignette vignette: > - %\VignetteEngine{R.rsp::md} + %\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{Multivariate RMSE} %\usepackage[utf8]{inputenc} --- @@ -18,7 +18,6 @@ To run this vignette, the next R packages should be installed and loaded: ```r library(s2dverification) library(RColorBrewer) -library(zeallot) ``` @@ -68,14 +67,18 @@ obs <- "erainterim" Using the `CST_Load` function from **CSTool package**, the data available in our data store can be loaded. The following lines show how this function can be used. Here, the data is loaded from a previous saved `.RData` file: Ask nuria.perez at bsc.es for the data to run the recipe. -```r +```r +require(zeallot) + glosea5 <- list(path = '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc') + c(exp_T, obs_T) %<-% CST_Load(var = temp, exp = list(glosea5), obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4, latmin = 25, latmax = 75, lonmin = -20, lonmax = 70, output = 'lonlat', nprocs = 1, storefreq = "monthly", sampleperiod = 1, nmember = 9, method = "bilinear", grid = paste("r", grid, sep = "")) + c(exp_P, obs_P) %<-% CST_Load(var = precip, exp = list(glosea5), obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4, diff --git a/vignettes/RainFARM_vignette.md b/vignettes/RainFARM_vignette.Rmd similarity index 93% rename from vignettes/RainFARM_vignette.md rename to vignettes/RainFARM_vignette.Rmd index 7b9817d19e5344bd474029d2e2f409e08d9c922c..90a882c871247a80d4aea561e95f92b019254737 100644 --- a/vignettes/RainFARM_vignette.md +++ b/vignettes/RainFARM_vignette.Rmd @@ -2,18 +2,16 @@ title: "Rainfall Filtered Autoregressive Model (RainFARM) precipitation downscaling" author: "Jost von Hardenberg (ISAC-CNR)" date: "26/03/2019" -output: - html_document: default - pdf_document: default +output: rmarkdown::html_vignette vignette: > + %\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{RainFARM} %\usepackage[utf8]{inputenc} - %\VignetteEngine{R.rsp::md} --- @@ -32,32 +30,33 @@ An efficient implementation of RainFARM is provided for CSTools by the `CST_Rain In order to run the examples in this vignette, the *CSTools* package and some other support R packages need to be loaded by running: -```r +```{r} install.packages('CSTools') library(CSTools) ``` We use test data provided by CSTools to load a seasonal precipitation forecast: -```r -exp %<-% lonlat_prec +```{r} +library(CSTools) +exp <- lonlat_prec ``` This gives us a CSTools object `exp`, containing an element `exp$data` with dimensions: -```r +```{r} dim(exp$data) -dataset member sdate ftime lat lon - 1 6 3 31 4 4 +#dataset member sdate ftime lat lon +# 1 6 3 31 4 4 ``` There are 6 ensemble members available in the data set, 3 starting dates and 31 forecast times, which refer to daily values in the month of March following starting dates on November 1st in the years 2010, 2011, 2012. Please notice that RainFARM (in this version) only accepts square domains, possibly with an even number of pixels on each side, so we always need to select an appropriate cutout. Also, there are time and memory limitations when a large ensemble of downscaled realizations is generated with RainFARM, so that selecting a smaller target area is advised. On the other hand, if spectral slopes are to be determined from the large scales we will still need enough resolution to allow this estimation. In this example we have preselected a 4x4 pixel cutout at resolution 1 degree in a smaller area lon=[6,9], lat=[44,47] covering Northern Italy. -#```r -#ilon <- which ( exp$lon %in% 5:12 ) -#ilat <- which ( exp$lat %in% 40:47 ) -#exp$data <- exp$data[ , , , , ilon, ilat, drop=FALSE] -#exp$lon <- exp$lon[ilon] -#exp$lat <- exp$lat[ilat] -#``` +```{r} +ilon <- which ( exp$lon %in% 5:12 ) +ilat <- which ( exp$lat %in% 40:47 ) +exp$data <- exp$data[ , , , , ilon, ilat, drop=FALSE] +exp$lon <- exp$lon[ilon] +exp$lat <- exp$lat[ilat] +``` ### Standard downscaling without climatological weights @@ -65,24 +64,25 @@ Our goal is to downscale with RainFARM these data from the resolution of 1 degre RainFARM can compute automatically its only free parameter, i.e. the spatial spectral slope, from the large-scale field (here only with size 4x4 pixel, but in general we reccomend selecting at least 8x8 pixels). In this example we would like to compute this slope as an average over the _member_ and _ftime_ dimensions, while we will use different slopes for the remaining _dataset_ and _sdate_ dimensions (a different choice may be more appropriate in a real application). To obtain this we specify the parameter `time_dim = c("member", "ftime")`. The slope is computed starting from the wavenumber corresponding to the box, `kmin=1`. We create 3 stochastic realizations for each dataset, member, starting date and forecast time with `nens=5`. The command to donwscale and the resulting fields are: -```r +```{r} exp_down <- CST_RainFARM(exp, nf=20, kmin = 1, nens = 3, time_dim = c("member", "ftime")) dim(exp_down$data) - dataset member realization sdate ftime lat lon - 1 6 3 3 31 80 80 +# dataset member realization sdate ftime lat lon +# 1 6 3 3 31 80 80 str(exp_down$lon) - num [1:80] 5.53 5.58 5.62 5.67 5.72 ... +# num [1:80] 5.53 5.58 5.62 5.67 5.72 ... str(exp_down$lat) - num [1:80] 47.5 47.4 47.4 47.3 47.3 ... +# num [1:80] 47.5 47.4 47.4 47.3 47.3 ... ``` The function returns an array `exp_down$data` with the additional "realization" dimension for the stochastic ensemble with 3 members. The longitudes and latitudes have been correspondingly interpolated to the finer resolution. Alternatively we could have used the "reduced" function `RainFARM` which accepts directly a data array (with arbitrary dimensions, provided a longitude, a latitude and a "time" dimension exist) and two arrays to describe longitudes and latitudes: -```r -downscaled <- RainFARM(exp$data, exp$lon, exp$lat, nf = 20, kmin = 1, nens = 3, +```{r} +downscaled <- RainFARM(exp$data, exp$lon, exp$lat, + nf = 20, kmin = 1, nens = 3, time_dim = c("member", "ftime")) ``` The resulting array has the same dimensions as the `$data` element in the result of `CST_RainFARM`. @@ -93,7 +93,7 @@ Each instant and each realization will of course be different, but let's plot an png("Figures/RainFARM_fig1.png", width = 640, height = 365) par(mfrow=c(1,2)) --> -```r +```{r} a <- exp$data[1, 1, 1, 17, , ] * 86400 * 1000 a[a > 60] <- 60 image(exp$lon, rev(exp$lat), t(apply(a, 2, rev)), xlab = "lon", ylab = "lat", @@ -119,13 +119,13 @@ The area of interest in our example presents a complex orography, but the basic Suitable climatology files could be for example a fine-scale precipitation climatology from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local high-resolution gridded climatology from observations, or a reconstruction such as those which can be downloaded from the WORLDCLIM (http://www.worldclim.org) or CHELSA (http://chelsa-climate.org) websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://www.gdal.org). We will assume that a copy of the WORLDCLIM precipitation climatology at 30 arcseconds (about 1km resolution) is available in the local file `medscope.nc`. From this file we can derive suitable weights to be used with RainFARM using the `CST_RFWeights` functions as follows: -```r +```{r} ww <- CST_RFWeights("./worldclim.nc", nf = 20, lon = exp$lon, lat = exp$lat) ``` The result is a two-dimensional weights matrix with the same `lon`and `lat` dimensions as requested. The weights (varying around an average value of 1) encode how to distribute differently precipitation in each stochastic realization of RainFARM. We call again `CST_RainFARM()`, this time passing climatological weights: -```r +```{r} exp_down_weights <- CST_RainFARM(exp, nf = 20, kmin = 1, nens = 3, weights = ww, time_dim = c("member", "ftime")) ``` @@ -134,8 +134,8 @@ We plot again the same realization as before in Fig. 2 (left panel). The WORLDCL From a single realization and time it is not possible to see that a more realistic precipitation has been achieved, but this becomes clear if we compare the climatological average over all 3 stochastic ensemble realizations, over all forecast times and over all forecast times. The resulting plot in panels Fig2b,c show that the RainFARM climatology downscaled without weights presents on average a very coarse structure, comparable with that of the original fields, while when using the weights a much more realistic distribution is achieved. -