Minimal dowsncaling example
The aim of this example is to show the basic steps for loading observations (station data) and model data (NCEP reanalysis) and operate with them in a statistical modelling framework. We illustrate how to establish a statistical link (a linear model) of spring (MAM) precipitation using the reference period 1981-2000 between the observed precipitation at the Barcelona Airport weather station (NE Spain) and the precipitation field given by the NCEP/NCAR reanalysis for that particular grid cell. Both datasets are built-in in the downscaleR toolbox, and are located in the corresponding sub-folders of the inst directory.
For the sake of simplicity, we avoid performing cross-validation or any other form of validation to avoid distracting from the essence of the message.
First of we load observations and model precipitation data for the selected location knowing the (approximated) geographical coordinates in degrees of Barcelona (2E, 41N). The necessary information can be obtained using the inventory and loading functions for station and gridded datasets.
> pr.obs <- loadObservations(source.dir="inst//datasets//observations//GSN_Iberia", var="precip", lonLim=2, latLim=41, season=3:5, years=1981:2000)
> pr.ncep <- loadGridDataset(dataset = "inst//datasets//reanalysis//Iberia_NCEP//Iberia_NCEP.ncml", var="tp", lonLim = 2, latLim = 41, season=3:5, years=1981:2000)
Next we inspect the loaded time series
> par(mfrow = c(1,2))
> plot(pr.obs$time$Start, pr.obs$Data$SP000008181, ty = 'l', xlab = "time", ylab = "Precip (mm/d)")
> lines(pr.ncep$Dates$Start, pr.ncep$Data, col = "red")
> legend("topright", c("Observed", "Reanalysis"), lty = 1, col = c(1,2))
> plot(pr.obs$Data$SP000008181, pr.ncep$Data, xlab = "Observed", ylab = "Reanalysis")
Next the linear model is fit using both data series:
> lm1 <- lm(pr.obs$Data$SP000008181 ~ pr.ncep$Data)
> summary(lm1)
Call:
lm(formula = pr.obs$Data$SP000008181 ~ pr.ncep$Data)
Residuals:
Min 1Q Median 3Q Max
-16.89 -0.78 -0.78 -0.78 63.82
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.78020 0.12131 6.431 1.61e-10 ***
pr.ncep$Data 0.93476 0.05753 16.249 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.992 on 1838 degrees of freedom
Multiple R-squared: 0.1256, Adjusted R-squared: 0.1251
F-statistic: 264 on 1 and 1838 DF, p-value: < 2.2e-16
Model predictions are obtained with the predict method for lmobjects
> predicted <- predict(lm1, pr.ncep)
Observed and predicted series are next plotted:
> plot(pr.obs$time$Start, pr.obs$Data$SP000008181, ty = 'l', xlab = "time", ylab = "Precip (mm/d)")
> lines(pr.obs$time$Start, predicted, col = "green")
> # Extraction of r squared value
> r2 <- summary(lm1)$adj.r.squared
> mtext(bquote(italic(r)^2 == .(format(r2, digits = 2))), cex = 1.5)
> legend("topright", c("Observed", "Predicted"), lty = 1, col = c(1,3))
