Changes between Version 2 and Version 3 of udg/ecoms/RPackage/examples/verification


Ignore:
Timestamp:
May 12, 2016 4:19:53 PM (6 years ago)
Author:
juaco
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • udg/ecoms/RPackage/examples/verification

    v2 v3  
    55
    66We will illustrate the verification of these predictions data against the observational gridded datasets WATCH Forcing
    7 Dataset-ERA-Interim (WFDEI, `dataset = "WFDEI"`), also available via the ECOMS-UDG. To this aim, we will use the tools developed within the projects SPECS and EUPORIAS. In particular, we will use the verification routines available in the R package `SpecsVerification` (available on CRAN). However, instead of using them directly, we will use the user-friendly interface implemented in package `easyVerification`, via the wrapper function `veriApply`.
     7Dataset-ERA-Interim (WFDEI, `dataset = "WFDEI"`), also available via the ECOMS-UDG. To this aim, we will use the tools developed within the projects SPECS and EUPORIAS. In particular, we will use the verification routines available in the R package `SpecsVerification` (available on CRAN). However, instead of using them directly, we will use the user-friendly interface implemented in package `easyVerification`, via the wrapper function `veriApply` (To find out more on the functionality in the `easyVerification` package, please read the vignette with the R instruction `vignette("easyVerification")`)
    88
    9 {{{
     9== Package loading/install
     10
     11We first load (and install if necessary) the required libraries. `loadeR.ECOMS` and `downscaleR` are loaded for data loading and manipulation respectively:
     12
     13{{{#!text/R
     14library(loadeR.ECOMS)
     15library(downscaleR)
     16}}}
     17
     18The two packages related with verification are next installed if not already present:
     19
     20
     21{{{#!text/R
     22# Attempting the installation of SpecsVerification
     23if (!require("SpecsVerification")) {
     24      install.packages("SpecsVerification")
     25}
     26# Attempting the installation of easyVerification
     27if (!require("easyVerification")) {
     28      if (!require(devtools)) {
     29            install.packages("devtools")
     30      }
     31      devtools::install_github("MeteoSwiss/easyVerification", build_vignettes=TRUE)
     32}
    1033}}}
    1134
    1235
    13 {{{#!text
     36== Data loading from the ECOMS-UDG
     37
     38{{{#!text/R
    1439tx.forecast <- loadECOMS(dataset = "CFSv2_seasonal",
    1540                         var = "tasmax",
     
    3358}}}
    3459
     60[[Image(image-20160512-135657.png)]]
    3561
    36 [[Image(image-20160512-135657.png)]]
     62
     63{{{#!text/R
     64tx.obs <- loadECOMS(dataset = "WFDEI",
     65                    var = "tasmax",
     66                    lonLim = c(-10 ,15),
     67                    latLim = c(35, 50),
     68                    season = 6:8,
     69                    years = 1991:2000)
     70## [2016-05-12 14:03:40] Defining homogeneization parameters for variable "tasmax"
     71## [2016-05-12 14:03:40] Opening dataset...
     72## [2016-05-12 14:03:42] The dataset was successfuly opened
     73## [2016-05-12 14:03:42] Defining geo-location parameters
     74## [2016-05-12 14:03:42] Defining time selection parameters
     75## [2016-05-12 14:03:42] Retrieving data subset ...
     76## [2016-05-12 14:03:52] Done
     77
     78plotMeanGrid(tx.obs)
     79}}}
     80
     81
     82[[Image(image-20160512-140857.png)]]
     83
     84== Data preprocessing: interpolation and aggregation using downscaleR
     85
     86=== Data interpolation
     87
     88To be able to validate the forecasts, we first have to interpolate either the
     89forecasts or the observations to have the data on a common grid. We interpolate the observations to the grid of the forecasts using the nearest neighbour algorithm:
     90
     91{{{#!text/R
     92tx.obsintp <- interpGrid(tx.obs, new.coordinates = getGrid(tx.forecast), method = "nearest")
     93## [2016-05-12 14:12:10] Calculating nearest neighbors...
     94## [2016-05-12 14:12:11] Performing nearest interpolation... may take a while
     95## [2016-05-12 14:12:11] Interpolating member 1 out of 4
     96## [2016-05-12 14:12:11] Interpolating member 2 out of 4
     97## [2016-05-12 14:12:11] Interpolating member 3 out of 4
     98## [2016-05-12 14:12:11] Interpolating member 4 out of 4
     99## [2016-05-12 14:12:13] Done
     100}}}
     101
     102=== Temporal aggregation
     103
     104Now that the predictions are in the same grid than the observations, we can proceed with the verification. However, before we go on, we compute seasonal averages of the forecasts and observations for validation of
     105seasonal average daily maximum temperature. This is easily undertaken using downscaleR's function `aggregateGrid` (in Linux and MacOS it is possible to speed-up the aggregation by using the parallelization option):
     106
     107{{{#!text/R
     108mn.tx.forecast <- aggregateGrid(tx.forecast, aggr.y = list(FUN = "mean"), parallel = TRUE)
     109## Parallel computing enabled
     110## Number of workers: 3
     111## [2016-05-12 14:30:17] Performing annual aggregation in parallel...
     112## [2016-05-12 14:30:18] Done.
     113
     114mn.tx.obsintp <- aggregateGrid(tx.obsintp, aggr.y = list(FUN = "mean"), parallel = TRUE)
     115##  Parallel computing enabled
     116##  Number of workers: 3
     117##  [2016-05-12 14:30:28] Performing annual aggregation in parallel...
     118##  [2016-05-12 14:30:30] Done.
     119}}}
     120
     121
     122Now we are ready to compute validation scores on the 3-monthly mean daily maximum temperature forecasts. We next compute several typical verification measures:
     123
     124== Verification measures using SpecsVerification and easyVerification
     125
     126=== Mean bias
     127
     128{{{#!text/R
     129bias <- veriApply("EnsMe",
     130                  fcst = mn.tx.forecast$Data,
     131                  obs = mn.tx.obsintp$Data,
     132                  ensdim = 1,
     133                  tdim = 2)
     134
     135}}}
     136
     137The object `bias` is a matrix with longitudes in columns and latitudes in rows. To do the plotting, we can use the coordinates of the reference grid of the CFS forecast. For convenience, we will use the function `image.plot` from package `fields` (this should be already installed because it is a dependency of `downscaleR`). In order to add a coastline map, we do a trick and call an internal function of `downscaleR` using the triple `:::` notation:
     138
     139{{{#!text/R
     140fields::image.plot(tx.forecast$xyCoords$x,
     141                   tx.forecast$xyCoords$y,
     142                   t(bias),
     143                   asp = 1, xlab = "", ylab = "", main = "Mean tmax bias - JJA")
     144downscaleR:::draw.world.lines()
     145
     146}}}
     147
     148[[Image(image-20160512-161841.png)]]
     149
     150
     151
     152
     153