Convert rasters to data frame with time stamp
Thank you Dominik, Vijay and Robers for all the valuable inputs.
Robert, the function you devised is really convenient. However, it fails when the observed and modelled raster objects have different number of layers.
Please notice that the obs and model objects on the example from the qmap packages have different lengths:
library(qmap)
data(obsprecip)
data(modprecip)
dim(obsprecip);dim(modprecip) # notice the difference in the lengths
#Fit a quantile mapping function to observed and modeled data
qm.fit <- fitQmap(obsprecip,modprecip,
method="QUANT",qstep=0.01)
#Perform bias correction on modeled data
qm <- doQmap(modprecip, qm.fit, type="tricub")
Now, an error is displayed when trying to apply your function on rasters with different "lengths" (number of layers):
library(raster)
library(qmap)
#Create a rasterStack similar to my data - same dimensions and layer names
r <- raster(ncol=60, nrow=60)
obs <- stack(lapply(1:408, function(x) setValues(r, runif(ncell(r)))))
mod <- stack(lapply(1:800, function(x) setValues(r, runif(ncell(r)))))
#Define bias-correction function
f <- function(obs, mod, ...) {
obs <- t(obs)
qm.fit <- fitQmap(obs, t(mod), method="QUANT",qstep=1/15)
t( doQmap(mod, qm.fit, type="tricub") )
}
# Test function on rasters with different # of layers
x <- f(head(obs), head(mod))
Error in t(doQmap(mod, qm.fit, type = "tricub")) :
error in evaluating the argument 'x' in selecting a method for function 't': Error in doQmapQUANT.matrix(x, fobj, ...) :
'ncol(x)' and 'nrow(fobj$par$modq)' should be eaqual
Any way to circumvent this?
Greetings,
-- Thiago V. dos Santos
PhD student
Land and Atmospheric Science
University of Minnesota
On Saturday, October 17, 2015 1:59 AM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
Thiago,
Therefore, I need to load data as rasters and iterate through all individual gridcells to create a data frame containing: date1, cell1, cell2, cell3, ..., cell3600 date2, cell1, cell2, cell3, ..., cell3600 date3, cell1, cell2, cell3, ..., cell3600... date408, cell1, cell2, cell3, ..., cell3600 then apply the fit function and finally convert the data back to a raster.
That would be easy enough by using, e.g., as.data.frame, as Dominik
and Vijay suggested, but I do not think that you _need_ to do that.
The manual says that you can use a matrix (no need for data.frame).
The dates are row.names, but I do not think they are used by the
algorithm, so you do not need to supply them (it would be possible).
In most cases like this, you can wrap the functions you need into a
short new function that can be passed to a raster function (calc or
overlay). In this case perhaps something like:
f <- function(obs, mod, ...) {
obs <- t(obs)
qm.fit <- fitQmap(obs, t(mod), method="QUANT",qstep=0.01)
t( doQmap(obs, qm.fit, type="tricub") )
}
try it
x <- f(head(s), head(s * 10))
x[, 1:5]
model <- s * 10
x <- overlay(s, model, fun=f)
Robert
On Fri, Oct 16, 2015 at 1:39 PM, Thiago V. dos Santos
<thi_veloso at yahoo.com.br> wrote:
Dear list,
Generally speaking, I love R. However, one of the things I like least in R is the need to interchange between the various data formats required by different packages.
I am trying to apply a bias-correction function on some gridded climate data. The qmap package has functions to perform bias correction on climate data, but the problem I am grasping with is that it requires data to be organized as data.frames:
library(qmap)
data(obsprecip)
data(modprecip)
#Fit a quantile mapping function to observed and modeled data
qm.fit <- fitQmap(obsprecip,modprecip,
method="QUANT",qstep=0.01)
#Perform bias correction on modeled data
qm <- doQmap(modprecip, qm.fit, type="tricub")
And that's all. But notice that both observed and modeled data in this example are data frames for different locations (Moss, Geiranger and Barkestad):
head(obsprecip)
MOSS GEIRANGER BARKESTAD 1-1-1961 0.1 0 0 2-1-1961 0.2 0 0 3-1-1961 0.9 0 0 4-1-1961 10.6 0 0 5-1-1961 1.5 0 0 6-1-1961 1.2 0 2
head(modprecip)
MOSS GEIRANGER BARKESTAD
2-1-1961 2.283 0.0000 3.177000
3-1-1961 2.443 10.8600 1.719000
4-1-1961 3.099 12.7300 6.636000
5-1-1961 0.000 9.7720 9.676000
6-1-1961 0.140 0.6448 7.110000
7-1-1961 13.470 3.3570 0.001107
Now, let's back to my problem. I have monthly precip data to which I want to apply the same function above, but my data is gridded:
library(raster)
#Create a rasterStack similar to my data - same dimensions and layer names
r <- raster(ncol=60, nrow=60)
s <- stack(lapply(1:408, function(x) setValues(r, runif(ncell(r)))))
names(s) <- paste0('X', seq(as.Date("1980/1/1"), by = "month", length.out = 408))
s
Therefore, I need to load data as rasters and iterate through all individual gridcells to create a data frame containing:
date1, cell1, cell2, cell3, ..., cell3600
date2, cell1, cell2, cell3, ..., cell3600
date3, cell1, cell2, cell3, ..., cell3600...
date408, cell1, cell2, cell3, ..., cell3600
then apply the fit function and finally convert the data back to a raster.
Any ideas on how to efficiently convert rasters to data frames containing their time stamp and then back to a raster again??
Any hint is much appreciated.
Greetings,
-- Thiago V. dos Santos
PhD student
Land and Atmospheric Science
University of Minnesota
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo