functionrasterqmap

Bias correction of raster climate data (qmap) - Error in formula

I am trying to apply a bias-correction function on some gridded climate data (observed and modelled projections).

The qmap package has functions to perform bias correction on climate data, I am trying to run the following

library(raster)
library(qmap)

#Create a rasterStack with observed and modelled data
r <- raster(ncol=20, nrow=20)
obs <- stack(lapply(1:100, function(x) setValues(r, runif(ncell(r)))))  #observed data
mod <- stack(lapply(1:100, function(x) setValues(r, runif(ncell(r)))))*2  #modelled data (i want this unbiased)

#bias-correction function 
f <- function(obs, mod, ...) {
  obs <- t(obs)
  qm.fit <- fitQmap(obs, t(mod), method="QUANT",qstep=0.01)
  t(doQmap(mod, qm.fit, type="linear") )
}

x <- overlay(obs, mod, fun=f)

I got the following error

Error in (function (x, fun, filename = "", recycle = TRUE, forcefun = FALSE, : cannot use this formula, probably because it is not vectorized

Can anybody help?

Thanks a million

Solution

  • With some light matrix reshaping, you can use the functions on matrices instead of rasters, which should work:

    library(qmap)
library(raster)


r <- raster(ncol=20, nrow=20)
obs <- stack(lapply(1:100, function(x) setValues(r, runif(ncell(r)))))  #observed data
mod <- stack(lapply(1:100, function(x) setValues(r, runif(ncell(r)))))*2  #modelled data (i want this unbiased)


qm.fit <- fitQmap(t(obs[]), t(mod[]), method="QUANT",qstep=0.01)

bias_corrected <- doQmap(t(mod[]), qm.fit, type="linear")

bias_corrected_arr <- as.array(t(bias_corrected))

#reshape array
dim(bias_corrected_arr) <- c(20,20,100)

# convert to rasterbrick
mod_Bcorrected <-setValues(brick(mod,values=FALSE),bias_corrected_arr)


# check output

> mod_Bcorrected
# class       : RasterBrick 
# dimensions  : 20, 20, 400, 100  (nrow, ncol, ncell, nlayers)
# resolution  : 18, 9  (x, y)
# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +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  : 2.136529e-03, 1.316528e-02, 4.380183e-03, 4.305932e-03, 5.293415e-03, 5.505609e-04, 9.494019e-05, 8.768495e-05, 8.171266e-04, 6.158992e-04, 3.248452e-03, 8.445294e-04, 2.856641e-03, 1.135400e-03, 3.255120e-03, ... 
# max values  :    0.9981889,    0.9995128,    0.9931242,    0.9975992,    0.9998942,    0.9986233,    0.9989456,    0.9952701,    0.9989686,    0.9958882,    0.9968628,    0.9975774,    0.9984780,    0.9910168,    0.9983078, ...