Changes between Version 8 and Version 9 of udg/ecoms/RPackage/examples/bias


Ignore:
Timestamp:
May 13, 2016 2:27:10 PM (6 years ago)
Author:
juaco
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • udg/ecoms/RPackage/examples/bias

    v8 v9  
    1919                       aggr.d = "min",
    2020                       aggr.m = "mean")
    21 }}}
    2221
    23 Some information messages will appear on-screen indicating the steps:
    24 
    25 {{{
    26 [2015-05-15 14:18:08] Defining homogeneization parameters for variable "tasmin"
    27 [2015-05-15 14:18:09] Defining geo-location parameters
    28 [2015-05-15 14:18:09] Defining initialization time parameters
    29 NOTE: Daily aggregation will be computed from 6-hourly data
    30 NOTE: Daily data will be monthly aggregated
    31 [2015-05-15 14:18:12] Retrieving data subset ...
    32 [2015-05-15 14:24:31] Done
     22## [2016-05-13 13:41:56] Defining homogeneization parameters for variable "tasmin"
     23## [2016-05-13 13:41:56] Opening dataset...
     24## [2016-05-13 13:41:56] The dataset was successfuly opened
     25## [2016-05-13 13:41:56] Defining geo-location parameters
     26## [2016-05-13 13:41:56] Defining initialization time parameters
     27## NOTE: Daily aggregation will be computed from 6-hourly data
     28## NOTE: Daily data will be monthly aggregated
     29## [2016-05-13 13:41:58] Retrieving data subset ...
     30## [2016-05-13 13:47:30] Done
    3331}}}
    3432
    3533
    36 Note the difference in size of the daily-aggregated data of the [https://meteo.unican.es/trac/wiki/udg/ecoms/RPackage/examples/continentalSelection in the previous example] (35.1 Mb) as compared to the new monthly-aggregated data size (1.2 Mb):
     34Note the difference in size of the daily-aggregated data of the [https://meteo.unican.es/trac/wiki/udg/ecoms/RPackage/examples/continentalSelection in the previous example] (35.1 Mb) as compared to the new monthly-aggregated data size (1.2 Mb). This is particularly important if our goal is to perform a validation of a large domain (possibly global), that can be easily handled on a monthly basis but that will cause memory problems using daily/sub-daily data.
    3735
    3836{{{#!text/R
    3937print(object.size(ex3.cfs), units = "Mb")
     38## 1.2 Mb
    4039}}}
    4140
    … …  
    5049                     years = 1991:2000,
    5150                     aggr.m = "mean")
    52 }}}
    53 {{{
    54 [2015-05-15 14:31:16] Defining homogeneization parameters for variable "tasmin"
    55 [2015-05-15 14:31:16] Defining geo-location parameters
    56 [2015-05-15 14:31:16] Defining time selection parameters
    57 NOTE: Daily data will be monthly aggregated
    58 [2015-05-15 14:31:17] Retrieving data subset ...
    59 [2015-05-15 14:31:39] Done
     51
     52## [2016-05-13 13:48:39] Defining homogeneization parameters for variable "tasmin"
     53## [2016-05-13 13:48:39] Opening dataset...
     54## [2016-05-13 13:48:39] The dataset was successfuly opened
     55## [2016-05-13 13:48:39] Defining geo-location parameters
     56## [2016-05-13 13:48:39] Defining time selection parameters
     57## NOTE: Daily data will be monthly aggregated
     58## [2016-05-13 13:48:39] Retrieving data subset ...
     59## [2016-05-13 13:49:08] Done
    6060}}}
    6161
    … …  
    7070[[Image(image-20150515-143320.png)]]
    7171
    72 Note that WFDEI provides data for land areas only, and its spatial resolution is finer than CFS (1º vs 0.5º). In order to compare both datasets, it is first necessary to put them in the same grid (i.e., to interpolate). We use bilinear interpolation to this aim, using the `downscaleR` function `interpGridData` in combination with the `getGrid` method, useful to recover the parameters defining the grid of a dataset to pass them to the interpolator:
     72Note that WFDEI provides data for land areas only, and its spatial resolution is finer than CFS (1º vs 0.5º). In order to compare both datasets, it is first necessary to put them in the same grid (i.e., to interpolate). We use bilinear interpolation to this aim, using the `downscaleR` function `interpGrid` in combination with the `getGrid` method, useful to recover the parameters defining the grid of a dataset to pass them to the interpolator:
    7373
    7474{{{
    7575#!text/R
    76 obs.regridded <- downscaleR::interpData(gridData = ex2.obs,
    77                                 new.coordinates = getGrid(ex2),
    78                                 method = "bilinear",
    79                                 parallel = TRUE)
     76obs.regridded <- downscaleR::interpGrid(grid = ex3.obs,
     77                                        new.coordinates = getGrid(ex3.cfs),
     78                                        method = "nearest")
     79## [2016-05-13 13:51:37] Calculating nearest neighbors...
     80## [2016-05-13 13:51:37] Performing nearest interpolation... may take a while
     81## [2016-05-13 13:51:37] Done
     82## Warning messages:
     83## 1: In downscaleR::interpGrid(grid = ex3.obs, new.coordinates = getGrid(ex2),  :
     84##   The new longitudes are outside the data extent
     85## 2: In downscaleR::interpGrid(grid = ex3.obs, new.coordinates = getGrid(ex2),  :
     86##   The new latitudes are outside the data extent
    8087}}}
    8188
    8289
    83 Note the use of the parallelization options to do the interpolation in parallel. Also, note the warnings reminding us that the extent of the input grid is wider than that from CFS. However, in this case we can safely ignore this warnings, since all the land areas we are interest in are within the CFS domain.
    84 
    85 {{{
    86 [2015-05-15 14:34:58] Performing bilinear interpolation... may take a while
    87 [2015-05-15 14:34:58] Done
    88 Warning messages:
    89 1: In interpGridData(gridData = ex3.obs, new.grid = getGrid(ex3.cfs),  :
    90   The new longitudes are outside the data extent
    91 2: In interpGridData(gridData = ex3.obs, new.grid = getGrid(ex3.cfs),  :
    92   The new latitudes are outside the data extent
    93 }}}
    94 
     90Note the warnings reminding us that the extent of the input grid is wider than that from CFS. However, in this case we can safely ignore this warnings, since all the land areas we are interest in are within the CFS domain.
    9591
    9692{{{
    … …  
    10096
    10197
    102 [[Image(image-20150515-143628.png)]]
     98[[Image(image-20160513-135405.png)]]
    10399
    104100
    105 After regridding, both model data and observations are in the same grid. We can compute the bias. First, we calculate the mean of WFDEI, which is the reference against which to compute the biases:
     101After regridding, both model data and observations are in the same grid. We can compute the bias.
    106102
    107 {{{
    108 #!text/R
    109 ref <- apply(obs.regridded$Data, MARGIN = c(3,2), mean, na.rm = TRUE)
    110 }}}
     103In order to compute the mean bias, we need to calculate the climatologies from both the observations, and the predictions. This can be easily accomplished using the function `climatology` in `downscaleR`. By default, it will compute the mean along the time dimension, but it is flexible and any other function could be defined by the user. In addition, the logical flag `by.member` allows to compute the climatology on each member by sepparate, rather than on the ensemble mean:
    111104
    112 The following lines of code compute the bias of each member w.r.t. the reference and plot them:
    113 
    114 {{{
    115 #!text/R
    116 # Now we compute the difference against each of the multimember spatial means:
    117 require(fields)
    118 n.members <- dim(ex3.cfs$Data)[1]
    119 par(mfrow = c(1,2))
    120 for (i in 1:n.members) {
    121       member <- apply(ex2$Data[i, , , ], MARGIN = c(3,2), mean, na.rm = TRUE)
    122       bias <- member - ref     
    123       image.plot(ex2$xyCoords$x, ex2$xyCoords$y, bias, xlab = "lon", ylab = "lat", asp = 1)
    124       title(paste("Bias member", i))
    125       world(add = TRUE)
    126 }
    127 par(mfrow = c(1,1)) # To reset the graphical window
     105{{{#!text/R
     106obs.clim <- climatology(obs.regridded)
     107## [2016-05-13 14:19:15] - Computing climatology...
     108## [2016-05-13 14:19:16] - Done.
     109pred.clim <- climatology(ex3.cfs, by.member = TRUE)
     110## [2016-05-13 14:19:31] - Computing climatology...
     111## [2016-05-13 14:19:31] - Done.
    128112}}}
    129113
    130114
    131 [[Image(image-20150515-144801.png)]]
     115Now, the biases can be computed, member by member, by just subtracting the observed climatology from the member climatologies. In this example, we create a new data array of the same dimensions than the predictions, include it in a new multimember grid object of the same characteristics as `ex3.cfs` and plot it:
    132116
     117{{{#!text/R
     118n.members <- dim(pred.clim$Data)[1]
     119arr <- pred.clim$Data
     120for (i in 1:n.members) {
     121      arr[i,,,] <- pred.clim$Data[i,,,] - obs.clim$Data[1,,,]
     122}
     123bias.grid <- ex3.cfs
     124bias.grid$Data <- arr
     125plotMeanGrid(bias.grid, multi.member = TRUE)
     126}}}
     127 
     128[[Image(image-20160513-142608.png)]]
     129
     130**NOTE**: a more elaborated example to compute bias as well as other verification metrics is provided is [../verification this EXAMPLE]
     131