santandermetgroup / downscaler

An R package for climate data bias correction and downscaling (part of the climate4R bundle)

What is downscaleR?

downscaleR is an R package for empirical-statistical downscaling focusing on daily data and covering the most popular approaches (bias correction, Model Output Statistics, Perfect Prognosis) and techniques (e.g. quantile mapping, regression, analogs, neural networks). This package has been conceived to work in the framework of both seasonal forecasting and climate change studies. Thus, it considers ensemble members as a basic dimension of the data structure. Find out more about this package at the downscaleR wiki.

This package is part of the climate4R bundle, formed by loadeR, transformeR, downscaleR and visualizeR. The recommended installation procedure is to use the install_github command from the devtools R package:

devtools::install_github(c("SantanderMetGroup/transformeR", "SantanderMetGroup/downscaleR"))

NOTE: Note that transformeR is a dependency for downscaleR. The utilities in transformeR were formerly part of downscaleR (up to v1.3-4). Since downscaleR v2.0-0, these are in transformeR and downscaleR is strictly aimed to statistical downscaling. Note that transformeR also includes illustrative datasets for the climate4rframework.

EXAMPLE: The following code trains three different downscaling methods (analogs, linear regression and neural networks) using principal components (explaining 95% of the variance for each variable) and visualizes the results (the illustrative station and reanalysis data for DJF included in the transformeR package is used in this example):

library(downscaleR)
data("VALUE_Iberia_tas") # illustrative datasets included in transformeR
y <- VALUE_Iberia_tas 
data("NCEP_Iberia_hus850", "NCEP_Iberia_psl", "NCEP_Iberia_ta850")
x <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
# calculating predictors
data <- prepareData(x = x, y = y,spatial.predictors = list(v.exp = 0.95)) 
# Fitting statistical downscaling methods (simple case, no cross-validation)
analog <- downscale.train(data, method = "analogs", n.analogs = 1)
regression <- downscale.train(data, method = "GLM",family = gaussian)
neuralnet <- downscale.train(data, method = "NN", hidden = c(10,5), output = "linear")
# Extracting the results for a particula station (Igueldo) for a single year (2000)
igueldo.2000 <- subsetGrid(y,station.id = "000234",years = 2000)
analog.2000 <- subsetGrid(analog$pred,station.id = "000234",years = 2000)
regression.2000 <- subsetGrid(regression$pred,station.id = "000234",years = 2000)
neuralnet.2000 <- subsetGrid(neuralnet$pred,station.id = "000234",years = 2000)
library(visualizeR)  # Data visualization utilities
temporalPlot(igueldo.2000, analog.2000, regression.2000, neuralnet.2000)

---

Reference and further information:

[General description of the downscaleR package] Bedia et al. (2020) Statistical downscaling with the downscaleR package (v3.1.0): contribution to the VALUE intercomparison experiment. Geosientific Model Development, 13, 1711–1735, https://doi.org/10.5194/gmd-13-1711-2020 Check out the companion notebooks GitHub.

[General description of the climate4R framework] Iturbide et al. (2019) The R-based climate4R open framework for reproducible climate data access and post-processing. Environmental Modelling and Software, 111, 42-54. https://doi.org/10.1016/j.envsoft.2018.09.009 Check out the companion notebooks for the two examples GitHub.

[Seasonal forecasting applications] Cofiño et al. (2018) The ECOMS User Data Gateway: Towards seasonal forecast data provision and research reproducibility in the era of Climate Services. Climate Services, 9, 33-43. http://doi.org/10.1016/j.cliser.2017.07.001

downscaler's Issues

Further explaination of the window argument in biasCorrection

I am having some difficulties in deeply understanding the machinery behind the optional "window" argument in the biasCorrection function.

A moving window is specified with two numbers: the window length and the window step. According to the wiki, the purpose is to obtain correction factors for each day of the year (doy). The window lenght is thus the number of days around each doy that is used to estimate the bias correction. The role of the window step is less clear: from my understanding if the step is 1, then the bias correction is done actually for each doy, it the step is n > 1, then the bias correction is done every n doy's. In practice, this will somewhat contradict the wiki page, since the correction will not be done for each doy, but actually at groups of doy's each consisting of n (step) days.

Is that correct? Thanks!

Potential bug when aggregate daily data during DJF

I have a SYS4 forecast of temperature for Winter DJF (subsetted for years):

> str(fcst)
List of 6
 $ Variable           :List of 2
  ..$ varName: chr "tas"
  ..$ level  : NULL
  ..- attr(*, "use_dictionary")= logi TRUE
  ..- attr(*, "description")= chr "2 metre temperature @ Ground or water surface"
  ..- attr(*, "units")= chr "degrees Celsius"
  ..- attr(*, "longname")= chr "2-meter air temperature"
  ..- attr(*, "daily_agg_cellfun")= chr "none"
  ..- attr(*, "monthly_agg_cellfun")= chr "none"
  ..- attr(*, "verification_time")= chr "none"
 $ Data               : num [1:15, 1:9388, 1:35, 1:50] 18.7 17.1 18.4 19.6 19.1 ...
  ..- attr(*, "dimensions")= chr [1:4] "member" "time" "lat" "lon"
 $ xyCoords           :List of 2
  ..$ x: num [1:50] -12 -11.25 -10.5 -9.75 -9 ...
  ..$ y: num [1:35] 32.2 33 33.7 34.5 35.2 ...
  ..- attr(*, "projection")= chr "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
 $ Dates              :List of 2
  ..$ start: chr [1:9388] "1983-12-01 00:00:00 GMT" "1983-12-01 06:00:00 GMT" "1983-12-01 12:00:00 GMT" "1983-12-01 18:00:00 GMT" ...
  ..$ end  : chr [1:9388] "1983-12-01 00:00:00 GMT" "1983-12-01 06:00:00 GMT" "1983-12-01 12:00:00 GMT" "1983-12-01 18:00:00 GMT" ...
  ..- attr(*, "subset")= chr "subsetYears"
  ..- attr(*, "season")= num [1:3] 12 1 2
 $ InitializationDates: chr [1:26] "1983-11-01 00:00:00 GMT" "1984-11-01 00:00:00 GMT" "1985-11-01 00:00:00 GMT" "1986-11-01 00:00:00 GMT" ...
 $ Members            : chr [1:15] "Member_1" "Member_2" "Member_3" "Member_4" ...
 - attr(*, "dataset")= chr "System4_seasonal_15"
 - attr(*, "source")= chr "ECOMS User Data Gateway"
 - attr(*, "URL")= chr "<http://meteo.unican.es/trac/wiki/udg/ecoms>"

When I perform the following operation:

  fcst_tas = aggregateGrid(fcst, aggr.d = list(FUN = mean))

the order of the time vector is broken:

> fcst_tas$Dates$start[1:100]
  [1] "1984-01-01 00:00:00 GMT" "1984-01-02 00:00:00 GMT" "1984-01-03 00:00:00 GMT"
  [4] "1984-01-04 00:00:00 GMT" "1984-01-05 00:00:00 GMT" "1984-01-06 00:00:00 GMT"
  [7] "1984-01-07 00:00:00 GMT" "1984-01-08 00:00:00 GMT" "1984-01-09 00:00:00 GMT"
 [10] "1984-01-10 00:00:00 GMT" "1984-01-11 00:00:00 GMT" "1984-01-12 00:00:00 GMT"
 [13] "1984-01-13 00:00:00 GMT" "1984-01-14 00:00:00 GMT" "1984-01-15 00:00:00 GMT"
 [16] "1984-01-16 00:00:00 GMT" "1984-01-17 00:00:00 GMT" "1984-01-18 00:00:00 GMT"
 [19] "1984-01-19 00:00:00 GMT" "1984-01-20 00:00:00 GMT" "1984-01-21 00:00:00 GMT"
 [22] "1984-01-22 00:00:00 GMT" "1984-01-23 00:00:00 GMT" "1984-01-24 00:00:00 GMT"
 [25] "1984-01-25 00:00:00 GMT" "1984-01-26 00:00:00 GMT" "1984-01-27 00:00:00 GMT"
 [28] "1984-01-28 00:00:00 GMT" "1984-01-29 00:00:00 GMT" "1984-01-30 00:00:00 GMT"
 [31] "1984-01-31 00:00:00 GMT" "1984-02-01 00:00:00 GMT" "1984-02-02 00:00:00 GMT"
 [34] "1984-02-03 00:00:00 GMT" "1984-02-04 00:00:00 GMT" "1984-02-05 00:00:00 GMT"
 [37] "1984-02-06 00:00:00 GMT" "1984-02-07 00:00:00 GMT" "1984-02-08 00:00:00 GMT"
 [40] "1984-02-09 00:00:00 GMT" "1984-02-10 00:00:00 GMT" "1984-02-11 00:00:00 GMT"
 [43] "1984-02-12 00:00:00 GMT" "1984-02-13 00:00:00 GMT" "1984-02-14 00:00:00 GMT"
 [46] "1984-02-15 00:00:00 GMT" "1984-02-16 00:00:00 GMT" "1984-02-17 00:00:00 GMT"
 [49] "1984-02-18 00:00:00 GMT" "1984-02-19 00:00:00 GMT" "1984-02-20 00:00:00 GMT"
 [52] "1984-02-21 00:00:00 GMT" "1984-02-22 00:00:00 GMT" "1984-02-23 00:00:00 GMT"
 [55] "1984-02-24 00:00:00 GMT" "1984-02-25 00:00:00 GMT" "1984-02-26 00:00:00 GMT"
 [58] "1984-02-27 00:00:00 GMT" "1984-02-28 00:00:00 GMT" "1984-02-29 00:00:00 GMT"
 [61] "1983-12-01 00:00:00 GMT" "1983-12-02 00:00:00 GMT" "1983-12-03 00:00:00 GMT"
 [64] "1983-12-04 00:00:00 GMT" "1983-12-05 00:00:00 GMT" "1983-12-06 00:00:00 GMT"
 [67] "1983-12-07 00:00:00 GMT" "1983-12-08 00:00:00 GMT" "1983-12-09 00:00:00 GMT"
 [70] "1983-12-10 00:00:00 GMT" "1983-12-11 00:00:00 GMT" "1983-12-12 00:00:00 GMT"
 [73] "1983-12-13 00:00:00 GMT" "1983-12-14 00:00:00 GMT" "1983-12-15 00:00:00 GMT"
 [76] "1983-12-16 00:00:00 GMT" "1983-12-17 00:00:00 GMT" "1983-12-18 00:00:00 GMT"
 [79] "1983-12-19 00:00:00 GMT" "1983-12-20 00:00:00 GMT" "1983-12-21 00:00:00 GMT"
 [82] "1983-12-22 00:00:00 GMT" "1983-12-23 00:00:00 GMT" "1983-12-24 00:00:00 GMT"
 [85] "1983-12-25 00:00:00 GMT" "1983-12-26 00:00:00 GMT" "1983-12-27 00:00:00 GMT"
 [88] "1983-12-28 00:00:00 GMT" "1983-12-29 00:00:00 GMT" "1983-12-30 00:00:00 GMT"
 [91] "1983-12-31 00:00:00 GMT" "1985-01-01 00:00:00 GMT" "1985-01-02 00:00:00 GMT"
 [94] "1985-01-03 00:00:00 GMT" "1985-01-04 00:00:00 GMT" "1985-01-05 00:00:00 GMT"
 [97] "1985-01-06 00:00:00 GMT" "1985-01-07 00:00:00 GMT" "1985-01-08 00:00:00 GMT"
[100] "1985-01-09 00:00:00 GMT"

Minor issue with subsetYears

In the documentation there is:

Logical (default is TRUE). Drop or keep dimensions of length 1

but there is no default.

> d = subsetYears(d, years = 1980:2005)
Error in asub.default(grid$Data, time.ind, grep("time", dimNames), drop = drop) : 
  argument "drop" is missing, with no default

Furthermore, there is a warning when executing the function:

> d = subsetYears(d, years = 1980:2005, drop = T)
Warning message:
In if (yrs == unique(yrs)) { :
  the condition has length > 1 and only the first element will be used

I'm running the last master release.

subsetSpatial with subsetting coordinates outside the extent

An exception in subsetSpatial for subsetting coordinates outside the grid extent is needed (also the corresponding warning message)

Bug in detrendGrid

Problem description:

tas.obs<-detrendGrid(temporal.tas.obs)
[2016-05-13 09:24:55] - Detrending...
[2016-05-13 09:26:20] - Done.
Error in aperm.default(arr, perm = match(newdim, refdim)) :
  value out of range in 'perm'
In addition: Warning message:
In newdim != refdim :
  longer object length is not a multiple of shorter object length

Reproducibility
get file at /oceano/gmeteo/users/maru/projects_R/specs/temporal.tas.obs_E-OBS.rda

plotEOF doesn't work with RStudio

If I try the first example into plotEOF documentation page:

# Winter temperature at 850 mb isobaric surface pressure level is loaded (period 1981--2010):
data(iberia_ncep_ta850)

PCA analysis, retaining the PCs that explain 95% of the total variance:

pca <- prinComp(iberia_ncep_ta850, v.exp = .95)

Plot of all EOFs

plotEOF(pca)

I get an error Error in plot.new() : figure margins too large

If I run the same code from terminal it works smoothly.

R.version
_
platform x86_64-apple-darwin13.1.0
arch x86_64
os darwin13.1.0
system x86_64, darwin13.1.0
status
major 3
minor 1.0
year 2014
month 04
day 10
svn rev 65387
language R
version.string R version 3.1.0 (2014-04-10)
nickname Spring Dance

Bias Correction for Vector data

Hello,
I have used isimip bias correction by following the procedure given in example. But it gives the results in form of same grids for whole area defined by latlon limits. I would really appreciate if someone can guide me if this method can be used in such a way that it returns the corrected values only for station grids instead of whole area.
For example I have 50 stations for temperature/ precipitation and I have 50 model times series data already extracted for the grids in which station falls, the bias correction should return the 50 times series.
Thank you in advance

Data has wrong time dimension when using delta biasCorrection

Here, I have an observational dataset (obs) and a forecast (fcst).

> str(obs)
List of 4
 $ Variable:List of 2
  ..$ varName: chr "var167"
  ..$ level  : NULL
  ..- attr(*, "is_standard")= logi FALSE
  ..- attr(*, "units")= chr "undefined"
  ..- attr(*, "longname")= chr "undefined"
  ..- attr(*, "daily_agg_cellfun")= chr "none"
  ..- attr(*, "monthly_agg_cellfun")= chr "none"
  ..- attr(*, "verification_time")= chr "none"
 $ Data    : num [1:72, 1:35, 1:50] 292 294 295 294 295 ...
  ..- attr(*, "dimensions")= chr [1:3] "time" "lat" "lon"
 $ xyCoords:List of 2
  ..$ x: num [1:50] -12 -11.25 -10.5 -9.75 -9 ...
  ..$ y: num [1:35] 32.2 33 33.7 34.5 35.2 ...
  ..- attr(*, "projection")= chr "LatLonProjection"
  ..- attr(*, "resX")= num 0.75
  ..- attr(*, "resY")= num 0.75
  ..- attr(*, "interpolation")= chr "nearest"
 $ Dates   :List of 2
  ..$ start: chr [1:72] "1984-06-30 18:00:00 GMT" "1984-07-31 18:00:00 GMT" "1984-08-31 18:00:00 GMT" "1985-06-30 18:00:00 GMT" ...
  ..$ end  : chr [1:72] "1984-06-30 18:00:00 GMT" "1984-07-31 18:00:00 GMT" "1984-08-31 18:00:00 GMT" "1985-06-30 18:00:00 GMT" ...
 - attr(*, "dataset")= chr "/opt/data/ERAIN-T2M/ERAIN-t2m-1983-2012.mon.EUROPE.nc"

> str(fcst)
List of 6
 $ Variable           :List of 2
  ..$ varName: chr "tas"
  ..$ level  : NULL
  ..- attr(*, "use_dictionary")= logi TRUE
  ..- attr(*, "description")= chr "2 metre temperature @ Ground or water surface"
  ..- attr(*, "units")= chr "degrees Celsius"
  ..- attr(*, "longname")= chr "2-meter air temperature"
  ..- attr(*, "daily_agg_cellfun")= chr "none"
  ..- attr(*, "monthly_agg_cellfun")=function (x, ...)  
  ..- attr(*, "verification_time")= chr "none"
 $ Data               : num [1:72, 1:35, 1:50] 21.8 23.3 23.2 19.1 21.1 ...
  ..- attr(*, "dimensions")= chr [1:3] "time" "lat" "lon"
 $ xyCoords           :List of 2
  ..$ x: num [1:50] -12 -11.25 -10.5 -9.75 -9 ...
  ..$ y: num [1:35] 32.2 33 33.7 34.5 35.2 ...
  ..- attr(*, "projection")= chr "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
 $ Dates              :List of 2
  ..$ start: chr [1:72(1d)] "1984-06-01 00:00:00 GMT" "1984-07-01 00:00:00 GMT" "1984-08-01 00:00:00 GMT" "1985-06-01 00:00:00 GMT" ...
  ..$ end  : chr [1:72(1d)] "1984-06-30 18:00:00 GMT" "1984-07-31 18:00:00 GMT" "1984-08-31 18:00:00 GMT" "1985-06-30 18:00:00 GMT" ...
 $ InitializationDates: chr [1:26] "1984-05-01 00:00:00 GMT" "1985-05-01 00:00:00 GMT" "1986-05-01 00:00:00 GMT" "1987-05-01 00:00:00 GMT" ...
 $ Members            : chr "Member_1"
 - attr(*, "dataset")= chr "System4_seasonal_15"
 - attr(*, "source")= chr "ECOMS User Data Gateway"
 - attr(*, "URL")= chr "<http://meteo.unican.es/trac/wiki/udg/ecoms>"

When I apply a delta bias correction with the command:

cal <- biasCorrection(y = obs,
                      x = fcst,
                      newdata = fcst,
                      cross.val = 'loocv',
                      method = "delta")

I got a grid with a wrong time dimension:

> str(cal)
List of 6
 $ Variable           :List of 2
  ..$ varName: chr "var167"
  ..$ level  : NULL
  ..- attr(*, "is_standard")= logi FALSE
  ..- attr(*, "units")= chr "undefined"
  ..- attr(*, "longname")= chr "undefined"
  ..- attr(*, "daily_agg_cellfun")= chr "none"
  ..- attr(*, "monthly_agg_cellfun")= chr "none"
  ..- attr(*, "verification_time")= chr "none"
 $ Data               : num [1:1656, 1:35, 1:50] 295 296 297 294 296 ...
 $ xyCoords           :List of 2
  ..$ x: num [1:50] -12 -11.25 -10.5 -9.75 -9 ...
  ..$ y: num [1:35] 32.2 33 33.7 34.5 35.2 ...
  ..- attr(*, "projection")= chr "LatLonProjection"
  ..- attr(*, "resX")= num 0.75
  ..- attr(*, "resY")= num 0.75
  ..- attr(*, "interpolation")= chr "nearest"
 $ Dates              :List of 2
  ..$ start: chr [1:72(1d)] "1984-06-01 00:00:00 GMT" "1984-07-01 00:00:00 GMT" "1984-08-01 00:00:00 GMT" "1985-06-01 00:00:00 GMT" ...
  ..$ end  : chr [1:72(1d)] "1984-06-30 18:00:00 GMT" "1984-07-31 18:00:00 GMT" "1984-08-31 18:00:00 GMT" "1985-06-30 18:00:00 GMT" ...
 $ InitializationDates: chr [1:26] "1984-05-01 00:00:00 GMT" "1985-05-01 00:00:00 GMT" "1986-05-01 00:00:00 GMT" "1987-05-01 00:00:00 GMT" ...
 $ Members            : chr "Member_1"
 - attr(*, "dataset")= chr "/opt/data/ERAIN-T2M/ERAIN-t2m-1983-2012.mon.EUROPE.nc"

This is not happening with other bias correction methods! It seems a bug due to CV, in fact without it, it works fine.

Error in if (Pth <= threshold) { : missing value where TRUE/FALSE needed_

Dear all,
I am using transformeR version 1.1.1 (2017-09-04) and downscaleR version 2.0.4 (2017-09-04).
I am getting the following error while performing bias correction on my files (y= observation data, x=model data) for training period

cal <- biasCorrection (y, x, x, precipitation = TRUE, method = "eqm")
[2017-09-06 01:17:28] Bias-correcting member 1 out of 1...
Error in if (Pth <= threshold) { : missing value where TRUE/FALSE needed

Although, I have interpolated both the files on a common grid and other attributes like missing value, are also kept same. I am not able to understand this error. Please help.

Thanks
Anubhav

array3Dto2Dmat seems quite inefficent

Hi,
I am implementing a set of analysis of gridded data based on PCA/SVD.
Looking at your code, I was wondering why array3Dto2Dmat is so complex.
It doesn't seem to perform any particular check. I implemented a trivial function that just permute the dimensions in the "right" order (assuming that we have just time, lat and lon) and change the number of dimensions from 3 to 2.
Here it is:

getX2 <- function(X) {
  X = aperm(X, match(c("time", "lat", "lon"), attr(X, 'dimensions')))
  dim(X) = c(dim(X)[1], dim(X)[2]*dim(X)[3])
  return(X)
}

I apologise, it's not really elegant. It was just to have a benchmark with your function.
X is a gridded data with 30x86x87 as dimensions.

mmm = microbenchmark(array3Dto2Dmat(X), getX2(X))
 > print(mmm)
Unit: milliseconds
              expr        min         lq       mean     median        uq      max neval cld
 array3Dto2Dmat(X) 206.937172 218.871806 230.222587 225.217348 234.11453 428.8247   100   b
          getX2(X)   1.275979   1.458082   2.083807   1.615363   2.23094  13.7102   100  a 

@jbedia what do you think?

Delta method is not correcting newdata (RCM climate projection) - BiasCorrection1D

Hello,

I've applied the methods "Linear Scaling, Delta, EQM and Variance to bias correction of temperature) with the function "BiasCorrection1D" (I have a single station to correct)

But the method "Delta" is not correcting new data.

When I run this block of code, I get the the correction for the "p" again, and not for "s".
I know it because "p" has 96 observations and "s" has 36 observations.

Bias Correction

library(downscaleR)

Loading mxtp observed data

o <- read.csv("o.csv", header = FALSE)
p <- read.csv("p.csv", header = FALSE)
s <- read.csv("s.csv", header = FALSE)

Bias correction (Delta)

delta <- downscaleR:::delta(o, p , s)

May someone help me?

Thanks in advance

A minor bug in gpqm function?

Dear all,

I was checking the R source code for gpqm function in biasCorrection.R and found an oddity from line 712. The line states that if auxF2=1, the adjusted precipitation values are obtained using max(o[indpareto]). If I have understood correctly, wrong indpareto values are now used since the same variable contains indices for p at this point, whereas indices calculated for o shoud be used.

Cheers,

Olle Räty

Example on biasCorrection documentation does not work

I'm using the last version of downscaleR.
I was trying to use the bias correction capabilities on my data but I got an error message. Then I tried with the example on biasCorrection documentation and...

> data(VALUE_Igueldo_tp)
> data(NCEP_Iberia_tp)
> y <- VALUE_Igueldo_tp
> x <- NCEP_Iberia_tp

> eqm1 <- biasCorrection(y = y, x = x, newdata = x, method = "eqm", extrapolation = "none", window = NULL)
Error in biasCorrection(y = y, x = x, newdata = x, method = "eqm", extrapolation = "none",  : 
  y and x do not have the same time series length

How to apply bias correction to RCP simulation by correcting historical simulation

Hai All,
I have doubt in biascorrection of future (say RCP) simulation of model by using historical simulation of model and observation.
Let me first correct the historical simulation bias by passing y as observation and x as historical simulation by using the below line.
cal.t <- biasCorrection (y = y.t, x = x.t, newdata = x.t, method = "eqm")

Now how should I apply the same bias correction to the future simulation (say RCP) outputs from the bias correction results of historical outputs ?

Thanks in advance.
Arulalan.T

datainventory

Hi there,

I find the tool very useful as it serves my purpose of bias correction. As I downloaded V2.0, and started to work on it as per Wiki, it seems like some of the functions as per wiki has not been updated in latest version like dataInventory, loadstationdata, etc. I further realized that it requires loadeR, but everytime installing these packages say that the package is not available for R version 3.3.2, while in description file it works for R > 3.0.0. Could you please suggest on this?

I am very new to this tool, so look forward to your valuable guidance. Moreover, pdf version of wiki with step by step guide would be much appreciated by new users like me. Thanking you, Nitesh

Latitudinal correction

A function for correcting by actual area is needed, to be applied on regular grids

Biascorrection using simple station and model output data (time_series)

Dear all ,

I am new to this package (downscaleR) and I found it very helpful. I am interested in Biascorrection and I have station and model output data(simulated and projected). I would like to start with one station and applied "biasCorrection(x, y, new, method = "eqm", theta = 0.95, window = NULL)", but this is not working. "Error in grepl(xx$Variable$varName, c("pr", "tp", "precipitation", "precip")) :
invalid 'pattern' argument". I need your help!!!
Thanks a lot!!!
Solomon

EQM formulation

Hello,

I was applying the EQM method in BiasCorrection1D (temperature) for a control Run period (without climate change scenarios), i.e.:
o <- y ;
p <- x;
s <- x;

Since the formulation presented in Amengual et al. (2012) for f=g=1 is:

EQM(i) = o(i) + Delta_m + Delta'(i)

And

Where:
EQM(i): ith ranked value projected corrected;
o(i): ith ranked value observed in Control Period;
Delta_m: mean delta change;
Delta'(i): individual delta change.

And

Delta(i) = s(i) - p(i)

s(i): ith ranked value raw future simulated;
p(i) ith ranked value raw control simulated;

Delta_m = mean(Delta(i))

Delta'(i) = Delta(i) - Delta_m

Since p = s,

Delta_i = 0
Delta_m = 0

And

EQM = o;

I do not understand why my values projected by the function EQM in R, do not are equal to the observed (o) to the Control Run.

I think it may be due some linear function to this, but I do not understand how.

Can anyone help me?
Thanks in advance

Include season as attribute

Needed to preserve the season information after seasonal aggregations

• getSeason only useful for hourly to monthly datasets, but not for seasonal to annual grids.
• Add a special attribute in Dates so season is always traceable
• Update climatology and getSeason accordingly...

Bug in interpData

I'm using the last devel version from github.

> data(iberia_ncep_ta850)
> par(mfrow = c(2,1))
> plotMeanGrid(iberia_ncep_ta850)
> # Bilinear interpolation to domain centered in Spain using a 0.5 degree resolution
> # in both X and Y axes
> t <- interpData(iberia_ncep_ta850, new.coordinates = list(x = c(-10,5,.5),
+                                                           y = c(36,44,.5)),
+                 method = "bilinear")
[2016-05-06 11:45:23] Performing multi-member bilinear interpolation... may take a while
[2016-05-06 11:45:23] Interpolating member 1 out of 1
Error in z[cbind(lx1 + 1, ly1)] : subscript out of bounds

error in bias correction when using a masked file as observation data

@jbedia

I am applying eqm method to bias correct my model data (x). my observation data (y) is masked to contain India and so is the model data (x). both are masked on a common boundary and have exactly same number of missing and non-missing values and the total number of grid points. but still I am getting this error on applying bias correction command

cal <- biasCorrection (y, x, x, precipitation = TRUE, method = "eqm")
[2017-09-08 21:33:07] Bias-correcting member 1 out of 1...
Error in if (Pth <= threshold) { : missing value where TRUE/FALSE needed

Please help. I am not able to understand this error. This problem is only with quantile mapping methods and does not occur when I use other methods like- multiplicative scaling or ptr. also when I replace the missing value in my observation file with some random meaningful constant, then the bias correction process is completed without any error but this is not a solution.
any help will be highly appreciated

biasCorrection with qqmap produces inconsistent results

In case a simple qqmap bias correcction is run with in-sample validation, the results from the biasCorrection are systematically biased.

prd.test <- biasCorrection(obs, prd, prd, 'qqmap')
plot(sort(obs$Data), sort(prd.test$Data)) ## all data points should fall on the 1-1 line

This is due to the quantile type=8 setting in calibrateProj. Setting instead type=4 produces consistent results with the probabilities returned in the ecdf function (these range from 1/n to 1). Alternatively, the quantiles could be directly computed (not using ecdf) corresponding to the quantile type setting.

getGrid won't work with rotated grids

This affects also the spatial subsetting capabilities using RCMs

Can I apply the Bias Correction methods in monthly series?

Can I apply the methods to monthly data?

I looked up in the code and the formulation of the methods. And I think it is possible, but I would like to confirm.

In case daily data is necessary, is it generated a monthly scale factor? As example, in Linear scaling, we could have one factor to each month (12 correction factors).

Can anyone help?

Missing dimensions in bias corrected grid data

I don't know if it's related with the issue #31
However,
after bias correcting a gridded field with:

cal <- biasCorrection(y = obs,
                      x = fcst,
                      newdata = fcst,
                     cross.val = 'kfold',
                     folds = as.list(1984:2007),
                      method = "eqm")

(btw, please add an example on the use of k-fold CV because I had to read the source code to understand what I had to put into folds)
If I check the dimensions of cal:

> str(cal$Data)
 num [1:2208, 1:35, 1:50] 293 293 294 294 294 ...

As you can see, the dimensions attribute is missing, creating problems with other functions like aggregateGrid, in fact:

> c_mon = aggregateGrid(cal, aggr.m = list(FUN = mean))
[2016-09-20 09:38:46] Performing monthly aggregation...
Error in apply(grid$Data, MARGIN = mar, FUN = function(x) { : 
  dims [product 1] do not match the length of object [3864000]

If I add manually the dimensions, everything works.

possible bug in implementing "loci" in biascorrect.R

Hi
On going through the R code for implementing loci, i think "s" is wrongly used instead of "p" at the following places:

        gcmr <- rev(sort(s))     (SHOULD be p]
        mgcm <- mean(s[which(s > Pgcm)], na.rm = TRUE)  [SHOULD be p]

fail to install this package from github

l run the follwing code in R:
library('devtools')
install_github(c("SantanderMetGroup/downscaleR"))

but it fail to install and the error message is displayed as follows：

• installing source package 'downscaleR' ...** R
  ** byte-compile and prepare package for lazy loading
  Error : object 'bindGrid.member' is not exported by 'namespace:transformeR'
  ERROR: lazy loading failed for package 'downscaleR'
• removing 'C:/Users/user/Documents/R/win-library/3.5/downscaleR'
  In R CMD INSTALL
  Installation failed: Command failed (1)

How to get netcdf file for output of bias correction

Hello all,
I am new to using downscaleR for bias correction.
I am able to apply different bias correction methods on my gridded RCM data using gridded observational data for e.g. using 'eqm' method as under,
cal <- biasCorrection (y = y, x = x, newdata = x, precipitation = TRUE, method = "eqm")
it is working fine and I am able to visualize the bias corrected results from "cal".
But I also want to get a netcdf file output of bias corrected data. How is it possible? Please give suggestions.

Regards
Anubhav

typo in onAttach messsage

el mensaje inicial de downscaleR (y loadeR también lo tiene así, creo) aparecen los símbolos >) en orden inverso. Es decir, el paréntesis cierra la llamada de la función devtools::install_github(), mientras que el > cierra la parte del mensaje que se refiere a ejecutar un código devtools::install_github(). Ahora está devtools::install_github()

plotMeanGrid / multiPlot problem

• Minor bug: Eliminate coastline display for rotated coordinates when multi.member=TRUE

The output of biasCorrection cannot be plotted with plotclimatology

After performed a biasCorrection of a seasonal forecast with biasCorrection I've tried to plot the output by using plotMeanGrid or plotClimatology but...

> plotClimatology(climatology(cal), 'coastline')
Error in aperm.default(obj$Data, perm) : 
  'perm' is of wrong length 4 (!= 7)

Issues with interpData and GS5

I use the devel version of downscaleR.
This is a quasi-reproducible example:

`YEARS = c(1996, 2009)
LEAD_TIME = 1
MONTHS = 6:8
SEASON = 'MAM'
SOURCE = 'GS5'
TARGET = 'ERA20C'

LABEL = sprintf('%s-%s-%s-%d-%d', SOURCE, TARGET, SEASON, YEARS[1], YEARS[2])

LOAD GS5

load(sprintf('GS5_S12_TAS_%s.Rdata', SEASON))
fcst = subsetGrid(fcst, latLim = c(40, 60), lonLim = c(0,20))

LOAD OBS

FNAME = 'data/ERA20C_moda_2t_1899_2010.seas.nc'
VARNAME = 't2m'

obs = loadGridData(FNAME,
var = VARNAME, dictionary = F,
season = MONTHS,
latLim = c(40, 60), lonLim = c(0,20),
years = seq(YEARS[1], YEARS[2]))

fcst_m = interpData(fcst, new.Coordinates = getGrid(obs), method = 'bilinear')
[2016-03-15 12:13:34] Performing multi-member bilinear interpolation... may take a while
[2016-03-15 12:13:34] Interpolating member 1 out of 12
Error in abind(list(x = c(0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20), y = c(40, :
can only supply one list-valued argument for ...`

Feature request: parallelized interpolation

Especially when working with ensembles interpolation may takes a lot of time, it would be great to have a parallel version.

Bias correcting one station

Hi,
I have an extracted gridded data for one station. I want to bias correct the extracted station using an observed data.
I used the biasCorrection1D command but it did not work.
I am not trying to bias-correct the grids but only bias- correcting and extracted grid data.
Is there an option to do this in downscaleR since it seems the function biasCorrection1D does not seem to work anymore.
Thanks!

using downscaleR into a R Markdown (R Notebook)

Have you ever used climate4R packages into a R Markdown environment? I am using RStudio and I have a piece of downscaling code that if runs inside a R Markdown block gives me a repeated warning:

data length [2] is not a sub-multiple or multiple of the number of rows [3]

while if outside R Markdown, it runs smoothly. Any experience on that?

repeated variable names in multigrids

Function ppModelSetUp for downscaling should differentiate grids with the same variable name in a multigrid (e.g. grids of different levels). Otherwise variable names in the multigrid need to be reassigned before downscaling.

Error loading "downscaleR"

hello @jbedia...Thank you for your great work.

I am trying to load the library "downdcaleR", but I am encountering this error:
> library(downscaleR)

Error: package or namespace load failed for ‘downscaleR’:
 object ‘bindGrid.member’ is not exported by 'namespace:transformeR'

Y also tried to reinstall the package, but it didn't work.

Thank you in advance.

Error while using isimip bias correction

Hello Everyone,
I'm very new with all this stuff so i really need help.
I'm trying to use isimip bias correction. i was able to load grid and station data.
First i tried interpolation by following the codes/example given at
https://github.com/SantanderMetGroup/downscaleR/blob/devel/R/isimip.R
so when i run this code,
obs <- interpGrid(obs, new.coordinates = getGrid(prd), method = "nearest")

The following error came,

NOTE: argument 'bilin.method' ignored for nearest neighbour interpolation
Error in obj$xyCoords$x : $ operator is invalid for atomic vectors

Then i tried to run code for isimip bias correction without interpolation but it also didn't work and gave following message.

simBC <- isimip (obs, prd, sim, threshold = 1)
Error in aperm.default(auxObs, match(attr(pred$Data, "dimensions")[setdiff(1:length(dimPred), :
value out of range in 'perm'
In addition: Warning message:
In if ((any(grepl(obs$Variable$varName, c("tas", "mean temperature", :
the condition has length > 1 and only the first element will be used

I would really appreciate if someone help me out.
Thank you!
Rizwan

Pending important wiki updates

In the wiki:

• Update regridding information in wiki page (points 2 grids, grids 2 points)
• Update temporal Aggregationinformation, and move to new page with full aggregateGrid examples

Help with downscale.cv

Hi, following the suggestion of @jbedia I am trying to implement a leave-one-year-out cross-validation directly with downscale.cv instead of implementing manually the procedure.

Currently, I do something like:

data <- prepareData(
        x = nc_data %>% subsetYears(years = train_years),
        y = stn_c4r %>% subsetYears(years = train_years)
      )
model1 <- downscale.train(data, method = "analogs", n.analogs = 3)
...

and it works. Instead, I tried with:

m1 <- downscale.cv(x = nc_data, y = stn_c4r, 
                         method = "analogs", folds = 4)

And I got:

[1] "fold: 1 --> calculating..."
[2018-05-30 14:30:21] - Scaling ...
Error in checkSeason(grid, base) : Inconsistent seasons among input grids

If I change the number of folds I get the same error but in a different fold calculation...
Do you have an updated (I'm using the devel channel) example to be used here?

Thank you in advance!

Wiki updates pending

Update perfect-prog wiki section and move to github wiki pages

Questions : How to access ensemble input netcdf file ?

Hai,

I have the following three queries. Could anyone answer these?

Question 1:
I need to access ECMWF model forecast data with the ensemble and do bias correction with joint distribution. I downloaded ECMWF TIGGE dataset locally and I want to access it, but I am unable to load with esemble axis. I am getting the following error "NOTE: The grid does not contain an Ensemble Axis: 'member' argument was ignored". I kept ensemble axis as first, followed by time, latitude, and longitude. Is there any further rules needed to be set in the nc file to access ensemble input ?

Question 2:
I want to do bias correct the most recent seasonal forecast based on reforecast and observation (using past 10 years). If y data is seasonal forecast for 10 years, x data is observation for same 10 years and newdata need to be same time length (i.e. 10 years) or could it be just 1 year (i.e. one season) ?

Question 3:
Do the biasCorrection function look for time axis (datewise) of y data and x data to apply bias correction or does it just ignore time axis then apply bias correction on length of y, x data ?

Thanks in advance.
Arulalan.T

Check plotMeanGrid for plotting climatologies

Need update to be consistent with the new 'climatology' function
Multimember climatologies yield an error

Dates list in multigrids

The Dates component created by makeMultiGrid should inherit the attributes from the dates of each component (e.g. the "season" attribute)

Error

Hi Dears,

Some months ago I used downscaleR, to made bias correction for time series obs and wrf outputs, now been updated the new code and found a problems, please can you help me?.

Macu_bc_obs <-Macu_bc$obs
Macu_fcst <- Macu_bc$fcst

y <- Macu_bc$obs
x <- Macu_fcst
newdata <- x
cal.t <- biasCorrection (y, x, newdata = NULL, method = "scaling",

•                      scaling.type = "additive")
  

  Error in obs$Variable : $ operator is invalid for atomic vectors

quickDiagnostics error: Error in subsetSpatial(field, lonLim, latLim)

Hello,
I'm following the guidelines on https://github.com/SantanderMetGroup/downscaleR/wiki/Simple-application-of-the-methods and get the following error on running quickDiagnostics.

I've tried modifying the location coordinates, but they are within the field as seen in plotMeanField(pred) image output.

Cheers,
Richard
‘downscaleR’ version 0.9-1
Rstudio: 0.98.1103
System: LinuxMint 17.1

Richard

``

precipitation

data(VALUE_Igueldo_tp)
obs <- VALUE_Igueldo_tp

plotMeanField(pred)
points(obs$xyCoords[,1], obs$xyCoords[,2], cex = 2, pch = 17)

precipitation

data(NCEP_Iberia_tp)
pred <- NCEP_Iberia_tp

cal <- biasCorrection (obs = obs, pred = pred, sim = pred, method = "eqm")
Warning messages:
1: In densfun(x, parm[1], parm[2], ...) : NaNs produced
2: In densfun(x, parm[1], parm[2], ...) : NaNs produced
3: In densfun(x, parm[1], parm[2], ...) : NaNs produced
4: In densfun(x, parm[1], parm[2], ...) : NaNs produced
5: In densfun(x, parm[1], parm[2], ...) : NaNs produced
6: In densfun(x, parm[1], parm[2], ...) : NaNs produced
7: In densfun(x, parm[1], parm[2], ...) : NaNs produced
8: In densfun(x, parm[1], parm[2], ...) : NaNs produced
9: In densfun(x, parm[1], parm[2], ...) : NaNs produced
str(cal)
List of 4
$ Variable:List of 2
..$ varName: chr "tp"
..$ level : NULL
..- attr(, "is_standard")= logi TRUE
..- attr(, "units")= chr "mm"
..- attr(, "longname")= chr "total precipitation amount"
..- attr(, "daily_agg_cellfun")= chr "sum"
..- attr(, "monthly_agg_cellfun")= chr "none"
..- attr(, "verification_time")= chr "DD"
$ Data : num [1:903, 1:7, 1:9] 0 0 2.6 14.7 13.2 14.1 0 30.7 7.7 6.1 ...
..- attr(, "dimensions")= chr [1:3] "time" "lat" "lon"
..- attr(, "threshold")= num 1
..- attr(, "correction")= chr "eqm"
$ xyCoords:List of 2
..$ x: num [1:9] -9.38 -7.5 -5.62 -3.75 -1.88 ...
..$ y: num [1:7] 33.3 35.2 37.1 39 41 ...
..- attr(, "projection")= chr "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
$ Dates :List of 2
..$ start: chr [1:903] "1990-12-01 GMT" "1990-12-02 GMT" "1990-12-03 GMT" "1990-12-04 GMT" ...
..$ end : chr [1:903] "1990-12-02 GMT" "1990-12-03 GMT" "1990-12-04 GMT" "1990-12-05 GMT" ...

• attr(*, "dataset")= chr "NCEP"
• attr(*, "source")= chr "ECOMS User Data Gateway"
• attr(*, "URL")= chr "http://meteo.unican.es/trac/wiki/udg/ecoms"

x <- cal$xyCoords$x[findInterval(x = obs$xyCoords[1], vec = cal$xyCoords$x)]
y <- cal$xyCoords$y[findInterval(x = obs$xyCoords[2], vec = cal$xyCoords$y)]

loc <- c(x,y)

quickDiagnostics(obs = obs, sim = pred, downscaled = cal, location = loc)

Error in subsetSpatial(field, lonLim, latLim) :
Subset longitude boundaries outside the current field extent:
(-2.0392)
In addition: Warning message:
In if (type == "daily") { :
the condition has length > 1 and only the first element will be used
``

NOT a issue. How to create input data (format used in example ) using netcdf file?

Error when applying an analog method on a single station with the downscale.predict function

I am using this package to carry out a downscaling between a gridded field and a set of stations. Everything works when:

data = prepareData(x = nc_data, 
                     y = stn_c4r) 
  
  model1 <- downscale.train(data, method = "analogs", n.analogs = 3)

but if I try:

newdata = prepareNewData(newdata = nc_data %>% subsetYears(years = Y),
                           data)
m1test = downscale.predict(newdata, model1)

I got an error:

Error in array(unlist(p), dim = dim(p[[1]])) : 
  'dims' cannot be of length 0

this is because in line 63 of downscale.predict you have:

p <- array(unlist(p), dim = c(n, dim(p[[1]])))

but in this case, given that stn_c4r is a single station, for the analog case in p[[1]] there is a numeric vector, leading to a NULL dimensions vector. Differently, when using GLM there you can find a matrix and then everything is working.
P.S. I am using the devel version of downscaleR

Error when interpolating field

I was trying to interpolate the variables obs_m (available here) with the following line:

> obs_m = interpGrid(obs_m, new.coordinates = getGrid(fcst_m), method = 'bilinear')
[2016-07-21 11:03:33] Performing bilinear interpolation... may take a while
Error in interp.old(x, y, z, xo = xo, yo = yo, ncp = 0, extrap = extrap,  : 
  duplicate data points: need to set 'duplicate = ..' 

The target grid is here:

> dput(getGrid(fcst_m))
structure(list(x = c(-178.875, 180), y = c(-39.8128523737712, 
70.0926903516245)), .Names = c("x", "y"), projection = "LatLonProjection", resX = 1.121484375, resY = 1.12148512985098)

the grid data seems ok to me, is this a bug?

Unable to install downscaleR on Windows7(64 bit)

Hello,
I'm trying to install downscaleR on Windows 7 (64 bit) and getting the following error:

library(devtools)
devtools::install_github(c("SantanderMetGroup/downscaleR.java","SantanderMetGroup/downscaleR"))
Downloading github repo SantanderMetGroup/downscaleR.java@master
Installing downscaleR.java
"C:/PROGRA1/R/R-311.1/bin/x64/R" --no-site-file --no-environ --no-save --no-restore CMD INSTALL
"C:/Users/tutul/AppData/Local/Temp/RtmpoloQGb/devtools14c057be6fba/SantanderMetGroup-downscaleR.java-357b9ec" --library="C:/Users/tutul/Documents/R/win-library/3.1"
--install-tests

• installing source package 'downscaleR.java' ...
  ** R
  ** inst
  ** preparing package for lazy loading
  Warning: package 'rJava' was built under R version 3.1.3
  ** help
  *** installing help indices
  ** building package indices
  ** testing if installed package can be loaded
  *** arch - i386
  Warning: package 'rJava' was built under R version 3.1.3
  *** arch - x64
  Warning: package 'rJava' was built under R version 3.1.3
• DONE (downscaleR.java)
  Downloading github repo SantanderMetGroup/downscaleR@master
  Installing downscaleR
  Skipping 1 packages not available: downscaleR.java
  "C:/PROGRA1/R/R-311.1/bin/x64/R" --no-site-file --no-environ --no-save --no-restore CMD INSTALL
  "C:/Users/tutul/AppData/Local/Temp/RtmpoloQGb/devtools14c036d6444c/SantanderMetGroup-downscaleR-0da25a5" --library="C:/Users/tutul/Documents/R/win-library/3.1"
  --install-tests
• installing source package 'downscaleR' ...
  ** R
  Error in parse(outFile) :
  C:/Users/tutul/AppData/Local/Temp/RtmpoloQGb/devtools14c036d6444c/SantanderMetGroup-downscaleR-0da25a5/R/biasCorrection.R:79:1: unexpected input
  78: #'
  79: <<
  ^
  ERROR: unable to collate and parse R files for package 'downscaleR'
• removing 'C:/Users/tutul/Documents/R/win-library/3.1/downscaleR'
  Error: Command failed (1)

How to fix this problem and how can I install this package??

Note that I downscaleR.java have been installed succesfully.

Please help.