In miquelcaceres/meteoland: Landscape Meteorology Tools
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(meteoland)
Introduction
Purpose
Reliable meteorological data are a basic requirement for hydrological and ecological studies at the landscape scale. Given the large spatial variation of meteorology over complex terrains, meteorological records from a single weather station are often not representative of entire landscapes. Studies made on multiple sites over a landscape require different meteorological series for each site; and other studies may require meteorological data series for all grid cells of a landscape, in a continuous way. In these cases, spatial correlation between the meteorology series of different sites or cells must be taken into account. For example, the sequence of days with rain of contiguous cells will normally be the same or very similar, even if precipitation amounts may differ. Finally, studies addressing the impacts of climate change on forests and landscapes require downscaling coarse-scale predictions of global or regional climate models to the landscape scale. When downscaling predictions for several locations in a landscape, spatial correlation of predictions is also important.
With the aim to assist research of climatic impacts on forests, the R package meteoland provides utilities to estimate daily weather variables at any position over complex terrains:
- Spatial interpolation of daily weather records from meteorological stations. }
- Statistical correction of meteorological data series (e.g. from climate models).}
Using meteoland package is easy, but some ideas and concepts must be addresed
to make the experience easy for new users. This vignette is intended as a working example to explain all the steps needed to get comfortable with meteoland workflow. We assume that the reader is more or less familiar with S4 classes of package sp. You will find detailed information regarding the package structure, its functions, and how calculations are done in vignette "Package meteoland".
Installing and loading the package
First of all, before starting to work with the package, we must install and load
the library. meteoland stable version can be found at CRAN
(https://CRAN.R-project.org/package=meteoland), and it can be
installed and loaded as any other R package:
install.packages("meteoland")
library(meteoland)
Alternatively, you can install the development version located at GitHub using the devtools package:
library(devtools)
install_github("miquelcaceres/meteoland", build_vignettes = TRUE)
library(meteoland)
Interpolation of daily weather
Spatial interpolation is required when meteorology for the area and period of interest cannot be obtained from local sensors. The nearest weather station may not have data for the period of interest or it may be located too far away to be representative of the target area.
Preparing weather data for interpolation
data("exampleinterpolationdata")
st_names = row.names(exampleinterpolationdata@MinTemperature)
st_coords = exampleinterpolationdata@coords
st_data = data.frame(X_UTM = st_coords[,1], Y_UTM = st_coords[,2],
elevation = exampleinterpolationdata@elevation,
row.names = st_names)
tmin = exampleinterpolationdata@MinTemperature
tmax = exampleinterpolationdata@MaxTemperature
prec = exampleinterpolationdata@Precipitation
rhum = exampleinterpolationdata@RelativeHumidity
Before starting using the package, you need to have access to the elevation (in m) and daily weather data series corresponding a set of locations (normally weather stations). Elevation is needed because interpolation routines perform corrections for differences in elevation between the reference locations and the target point. The initial format of your data will be different depending on the format used by your data provider (the package has also tools to access weather data). For our example, we will assume you have data from a set of 38 stations in your study area. On one hand, you should have a data.frame with the coordinates and elevation of each location:
str(st_data)
head(st_data)
On the other, you should have at least three matrices of meteorological data (one for minimum temperature, one for maximum temperature and the last one for precipitation) with stations in rows and dates in columns. In our example we also add relative humidity (in percent), so that other derived variables can be calculated:
dim(tmax)
dim(tmin)
dim(prec)
dim(rhum)
tmax[1:6,1:6]
Units should be in degrees Celsius for temperature and mm for precipitation.
Building an interpolation data object
Meteoland stores weather series for reference locations and interpolation parameters in a single object of class MeteorologyInterpolationData. There are several ways of building such objects, but we will first illustrate how to do it from the data we just presented.
First we need to create an object of class SpatialPoints (see package sp) with the spatial coordinates of our stations and the coordinate system (here UTM 31N):
sp = SpatialPoints(st_data[,c("X_UTM", "Y_UTM")],
proj4string = CRS("+proj=utm +zone=31 +ellps=WGS84 +datum=WGS84 +units=m +towgs84=0,0,0"))
head(sp)
We can now build an object MeteorologyInterpolationData using:
interpolator <- MeteorologyInterpolationData(sp, elevation = st_data$elevation,
MinTemperature = tmin,
MaxTemperature = tmax,
Precipitation = prec,
RelativeHumidity = rhum)
class(interpolator)
The resulting object is now ready to be used to perform interpolation on a set of target locations (see next section). We can inspect the amount of data in our interpolation object using function interpolation.coverage. For example, with:
spatial_coverage <- interpolation.coverage(interpolator, type = 'spatial')
head(spatial_coverage)
we obtain the number of non-missing observations for each weather station and variable. Similarly, we can ask the number of non-missing observations for each date and variable using:
temporal_coverage <- interpolation.coverage(interpolator, type = 'temporal')
head(temporal_coverage)
Interpolation parameters are also stored in the same object (see next subsection):
names(interpolator@params)
Interpolation basics and parameters
Package meteoland implements, with a few modifications, the daily weather interpolation and estimation algorithms that un- derpin the U.S. DAYMET dataset (Thornton et al., 1997; Thornton and Running, 1999).This approach, similar to inverse distance weighting, interpolates weather variables using trun- cated Gaussian filters, which consist in defining spatial weights $W(r)$ at radial distance $r$ from a target point $p$ using:
\begin{equation}
W(r) = e^{-\alpha \cdot (r/R_p)^2} - e^{-\alpha}
\end{equation}
if $r \leq R_p$ and $W(r) = 0$ otherwise.
Here $r$ is the radial distance from $p$, $R_p$ is the truncation distance and $\alpha$ is the shape parameter. The spatial convolution of this filter with a set of weather station locations results, for each target point, in a vector of weights associated with observations. The following figure illustrates the Gaussian filter for $R_p = 500$ and either $\alpha = 3.0$ (continuous line) or $\alpha = 6.25$ (dashed line):
r = 0:1000
R_p = 500
gf1 = exp(-3.0*((r/R_p)^2.0)) - exp(-3.0)
gf2 = exp(-6.25*((r/R_p)^2.0)) - exp(-6.25)
gf1[r>R_p] = 0
gf2[r>R_p] = 0
plot(r, gf1, type="l", ylab = "W(r)", xlab ="r")
lines(r, gf2, lty=2)
legend("topright", legend = c("alpha = 3", "alpha = 6.25"), lty=c(1,2), bty="n")
$R_p$ is automatically adjusted so that it has lower values in data-rich regions and is increased in data-poor regions. The method, however, requires the user to specify $N$, the average number of observations to be included for each target point. $R_p$ is then varied as a smooth function of the local density in such a way that this average is achieved over the spatial domain. It is important that the initial value of $R_p$
In meteoland estimation of $Rp$ is done once for each target point, variable and day. Interpolation of temperature includes a correction for the effects of elevation. More specifically, a weighted least-squares regression is used to assess the relationship between temperature differences and elevation differences in weather station data and this relationship is applied to elevation differences between weather stations and the target point. Interpolation of relative humidity is done after transforming it to dew-point temperature. No correction for elevation is performed during interpolation, but elevation effects arise when back-transforming dew-point temperature to relative humidity. Interpolation of daily precipitation is complicated by the need to predict both precipitation occurrence and, conditioned on this, precipitation amount. Thornton et al. (1997) defined a binomial predictor of spatial precipitation occurrence as a function of the weighted occurrence at surrounding weather stations. Conditional on precipitation occurrence, the interpolation routine predicts precipitation amount, where weighted least-squares regression is also used to account for elevation effects. Interpolation of wind is performed in three different ways depending on the information available. If only wind speed data is available, the spatial interpolation with Gaussian weights is used on wind scalars as described above. If weather station data includes wind direction, a polar average is cal- culated using Gaussian weights. Finally, if static wind fields are also available the interpolation routine first finds, for each weather station the wind field that best matches the observed vector. Then, the wind vectors extracted from the selected wind fields are averaged as before. The details of how $R_p$ is estimated and, in general, how interpolation is done are given in a PDF vignette called 'Meteorology'.
Interpolation parameters $\alpha$ and $N$ can be different for each variable to be interpolated. The following table lists all the interpolation parameters (see also function defaultInterpolationParameters())
| Paremeter | Default value | Definition |
| ------------- | ---------- | ----------------------------------- |
| initial_Rp | 140000 | Initial truncation radius |
| iterations | 3 | Number of station density iterations |
| alpha_MinTemperature | 3.0 | Gaussian shape parameter for minimum temperature |
| alpha_MaxTemperature | 3.0 | Gaussian shape parameter for maximum temperature |
| alpha_DewTemperature | 3.0 | Gaussian shape parameter for dew-point temperature |
| alpha_PrecipitationEvent | 5.0 | Gaussian shape parameter for precipitation events |
| alpha_PrecipitationAmount | 5.0 | Gaussian shape parameter for the regression of precipitation amounts |
| alpha_Wind | 3.0 | Gaussian shape parameter for wind |
| N_MinTemperature | 30 | Average number of stations with non-zero weights for minimum temperature |
| N_MaxTemperature | 30 | Average number of stations with non-zero weights for maximum temperature |
| N_DewTemperature | 30 | Average number of stations with non-zero weights for dew-point temperature |
| N_PrecipitationEvent | 5 | Average number of stations with non-zero weights for precipitation events |
| N_PrecipitationAmount | 20 | Average number of stations with non-zero weights for the regression of precipitation amounts |
| N_Wind | 2 | Average number of stations with non-zero weights for wind |
| St_Precipitation | 5 | Number of days for the temporal smoothing of precipitation |
| St_TemperatureRange | 15 | Number of days for the temporal smoothing of temperature range |
| pop_crit | 0.50 | Critical precipitation occurrence parameter |
| f_max | 0.6 | Maximum value for precipitation regression extrapolations (0.6 equals to a maximum of 4 times extrapolation) |
| wind_height | 10 | Wind measurement height (in m) |
Parameter St_Precipitation controls the temporal smoothing that is applied to weather station data to calibrate regression, while parameters pop_crit and f_max are also particular to the estimation of precipitation (see details in Thornton et al. 1997). Parameter St_TemperatureRange is used for the estimation of solar radiation (see details in PDF vignette called 'Meteorology').
A parameter that is particularly important to understand is initial_Rp, which specifies the initial radius for the truncated spatial Gaussian kernel. By default its value is:
interpolator@params$initial_Rp
The value of initial_Rp must be set in relation to the units of the spatial coordinates of weather data. Our data was in meters, so the default radius is 140 km. In general, the initial radius should be large enough to include a reasonable number of stations (~20-40), but the kernel radius is adjusted for each interpolation target point.
Calibration and cross-validation of the interpolation data
Once we already have weather stations data in shape, we can start calibrating the model in order to obtain the optimal parameters for the meteorological variables we want to interpolate. Parameter calibration has to be done for each variable separately. For example for minimum temperature:
tmin_cal <- interpolation.calibration(interpolator, variable = "Tmin",
N_seq = 20,
alpha_seq = seq(5, 10, by = 1),
verbose = TRUE)
This function returns an interpolation.calibration class object which contains several items:
- Numeric matrix with the mean absolute error (MAE) values for each combination of parameters $N$ and
$\alpha$.
- Miminum value found for MAE.
- Value for the $N$ parameter corresponding to the minumun MAE.
- Value for the $\alpha$ parameter corresponding to the minimum MAE.
- Matrix with the observed values.
- Matrix with the predicted values for the optimum parameter combination.
The result of the calibration needs to be manually stored in the interpolation params:
interpolator@params$N_MinTemperature = tmin_cal$N
interpolator@params$alpha_MinTemperature = tmin_cal$alpha
We strongly recommend conducting calibration exercises at least once for each variable and each data set used as reference for interpolation, and more than once if periods differ in the number of stations available.
Before using the object for interpolations, we also need to assess its performance. This is done by cross-validation in function interpolation.cv:
cv <- interpolation.cv(interpolator, verbose=T)
Cross-validation is perfomed by leave-one-out, which means that each target station is first removed from the data set and then the remaining stations are used to obtain interpolation estimates. In this way one maximizes the amount of information for estimates, while keeping them independent of the observed values in each station. The results of cross-validation can be inspected using a specific summary function:
summary(cv)
Interpolation on a grid
data("examplegridtopography")
elev = as(examplegridtopography[1:20,1:20], "SpatialGridDataFrame")["elevation"]
The target for weather interpolation in meteoland can be a set of points, pixels or a whole grid. Again, the initial format of data can be very different. Here we assume you have a small grid of 400 (20x20) cells of 1ha in size, with elevation data in form of class SpatialGridDataFrame (see method read.asciigrid in package sp):
summary(elev)
Note that the coordinate reference system needs to be the same as that of interpolator, which in this case it is. Before performing the interpolation over this grid, we need to reshape this data in a class called SpatialGridTopography:
sgt = SpatialGridTopography(as(elev, "SpatialGrid"), elevation = elev$elevation,
proj4string = elev@proj4string)
sgt
As you can see in the result, meteoland has calculated for us slope and aspect (both in degrees) from elevation data ^[Objects of class SpatialGridTopography can be initialized with user input values for slope and aspect too, but meteoland has its own routines when this are missing]. Slope and aspect are important for radiation calculations, which also requires relative humidity data. We can display elevation over the grid using:
spplot(sgt, "elevation")
Before we call the interpolation routine, we need to define the dates (i.e. days) for which we want weather to be interpolated, for example:
dates = as.Date(c("2001-02-03", "2001-06-03"))
Of course, we need to be sure that the interpolator object has data corresponding to this dates. We can check if there is any missing date using:
sum(!(dates %in% interpolator@dates))
The name of interpolation functions depend on the target spatial structure. For grids we need to use function interpolationgrid:
ml <- interpolationgrid(interpolator, sgt, dates)
This function works processing each date at a time. Since calculations can take some time, the console output shows the progress. The output of the function is an object of class SpatialGridMeteorology:
ml
We can display interpolated grids in a map using function spplot:
spplot(ml, 2, "MinTemperature")
spplot(ml, 2, "MaxTemperature")
Objects of class SpatialGridMeteorology include a list of data frames, one per date. We can access the interpolated data for a given date using:
df_1 = ml@data[[1]]
head(df_1)
Some columns are missing (e.g. wind speed) because we did not include weather station data regarding these variables. If we wants to add grid coordinates to this data frame, we can use the spatial information stored in the ml object:
sgdf_1 = SpatialGridDataFrame(grid = ml@grid, data = ml@data[[1]], proj4string = ml@proj4string)
summary(sgdf_1)
Results can be retrieved in this way and saved using write.asciigrid for its use outside R. The package also provides function for reading/writing NetCDFs (see function writemeteorologygrid) from SpatialGridMeteorology objects.
Interpolation on a set of points
spt = as(examplegridtopography, "SpatialPointsTopography")[c(36,228, 530, 2304)]
cc = coordinates(spt)
points_df = data.frame(X_UTM = cc[,1], Y_UTM = cc[,2], elevation = spt@data$elevation)
If you want to interpolate on a set of target locations, the starting point will normally be a data.frame. In our example this points come from the grid, but we have reshaped them so the starting format is familiar:
points_df
Analogously with the grid, we need to transform this data into an object SpatialPointsTopography:
spt = SpatialPointsTopography(as.matrix(points_df[,c("X_UTM", "Y_UTM")]),
elevation = points_df$elevation,
proj4string = CRS("+proj=utm +zone=31 +ellps=WGS84 +datum=WGS84 +units=m +towgs84=0,0,0"))
spt
In this case we only have elevation (in m), but slope and aspect should also be included if possible. Let us assume you want to interpolate on this points for the whole time series available in object interpolator. Since we are dealing with points, the function to interpolate is called interpolationpoints:
mp = interpolationpoints(interpolator, spt)
This function works processing one point at a time. The output of the function is an object of class SpatialPointsMeteorology:
mp
And the time series for a given can be plotted using function meteoplot. For example, we show here the precipitation series of point #1:
meteoplot(mp, 1, "Precipitation", ylab="Precipitation (mm)", xlab="")
Objects of class SpatialPointsDataFrame include a list of data frames, one per point. We can access one of them using:
df_1 = mp@data[[1]]
head(df_1)
This data frame can now be written into a file for its analysis outside R. The package also provides its own functions to write/read point meteorology data in different formats. If we are interested in inspecting the interpolation result by date, instead of by point, we can use function extractpointdates, which returns objects of class SpatialGridDataFrame:
dt_4 = extractpointdates(mp, as.Date("2001-01-04"), verbose = FALSE)
dt_4
Statistical correction of daily weather
Correcting the biases of a meteorological data series containing biases using a more accurate meteorological series is necessary when the more accurate series does not cover the period of interest and the less accurate series does. The less accurate series may be at coarser scale, as with climate model predictions or climate reanalysis data. In this case one can speak of statistical correction and downscaling. However, one may also correct the predictions of climate models using reanalysis data estimated at the same spatial resolution.
In the following example we will correct the predictions of Regional Climate Model (CCLM4-8-17; driving global model CNRM-CERFACS-CNRM-CM5) on the same area of the interpolation example. RCM data includes 3 model cells. Meteorological data covers an historical (reference) period (2000-2003) and a future (projection) period (year 2023), the latter simulated under rcp4.5 scenario.
Preparing weather data to be downscaled/corrected
data("examplecorrectiondata")
coords = examplecorrectiondata@coords
pt_coords = data.frame(long = as.numeric(coords[,1]), lat = as.numeric(coords[,2]))
proj_dates = examplecorrectiondata@dates
ref_data = examplecorrectiondata@reference_data
proj_data = examplecorrectiondata@projection_data
One needs several data items to perform downscaling and statistical correction. First, we need a matrix or data frame with the central coordinates of the RCM cells (here in longitude/latitude format):
pt_coords
Second, we need the uncorrected meteorological data (here RCM outputs) for both a reference period (that will be matched with our more accurate series) and a projection period (that our more accurate series does not cover). Both should be arranged in lists of data frames, one per RCM cell:
length(ref_data)
length(proj_data)
head(ref_data[[1]])
head(proj_data[[1]])
Note that these data frames should follow the conventions of meteoland for variable names and units. Finally, we need to specify the projection period, in this case one year:
proj_dates = seq(as.Date("2023-01-01"), as.Date("2023-12-31"), by="day")
Building the uncorrected data object
Statistical correction needs an object of class MeteorologyUncorrectedData, analogous to the MeteorologyInterpolationData of the previous section. To build this object we need first to express coordinates in an object SpatialPoints:
sp = SpatialPoints(pt_coords, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
sp
Assuming the weather data to be corrected is in the proper format, the object is constructed using:
uncorrected = MeteorologyUncorrectedData(sp, ref_data, proj_data, proj_dates)
Target point weather meteorology taken as reference
Downscaling/statistical correction is done on a set of target spatial points. meteoland will first find the RCM cell to which each point is nearest (using the central coordinates of RCM cells). Once this matching is done, statistical relations are build between the meteorology of RCM cells and that of target points for the reference period (here 2000-2003) which are used to correct RCM data for the projection period (year 2023). For this process to be done, we need the coordinates and weather series for the reference period. In our case, we will employ a copy the object of class SpatialPointsMeteorology that resulted from interpolation in the previous section:
historical = mp
historical
We named it historical to remind that we are using the historical period as reference and want to correct the RCM series for the projection period There are several ways to shape data into objects of class SpatialPointsMeteorology. For example, one can use function readmeteorologypoints.
Tuning correction parameters
Objects of class MeteorologyUncorrectedData have a slot with parameters:
uncorrected@params
Importantly, varmethods is a named list that specifies which correction method has to be used for each variable. Since we do not have reference data for wind speed, we must turn the method for wind speed to "none":
uncorrected@params$varmethods$WindSpeed="none"
Conducting statistical correction
Once we have all the data objects that we need (the hard part), conducting statistical correction is straightforward:
projected = correctionpoints(uncorrected, historical)
As mentioned above, correction proceeds point by point. The result is again an object of class SpatialPointsMeteorology:
projected
The following code displays the minimum/maximum temperatures before and after correction:
#Plot predicted mean temperature for point 1
meteoplot(projected, 1, "MinTemperature", ylab="Temperature (Celsius)", ylim=c(-5,40), col="blue")
meteoplot(projected, 1, "MaxTemperature", add=TRUE, col="red")
#Add uncorrected mean temperature data (cell #3)
lines(uncorrected@dates,
uncorrected@projection_data[[3]]$MinTemperature,
col="blue", lty=3)
lines(uncorrected@dates,
uncorrected@projection_data[[3]]$MaxTemperature,
col="red", lty=3)
legend("topright", legend=c("corrected","uncorrected", "Maximum", "Minimum"),
col=c("black","black", "red","blue"), lty=c(1,3,1,1), bty="n")
miquelcaceres/meteoland documentation built on May 8, 2019, 11:57 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(meteoland)
Introduction
Purpose
Reliable meteorological data are a basic requirement for hydrological and ecological studies at the landscape scale. Given the large spatial variation of meteorology over complex terrains, meteorological records from a single weather station are often not representative of entire landscapes. Studies made on multiple sites over a landscape require different meteorological series for each site; and other studies may require meteorological data series for all grid cells of a landscape, in a continuous way. In these cases, spatial correlation between the meteorology series of different sites or cells must be taken into account. For example, the sequence of days with rain of contiguous cells will normally be the same or very similar, even if precipitation amounts may differ. Finally, studies addressing the impacts of climate change on forests and landscapes require downscaling coarse-scale predictions of global or regional climate models to the landscape scale. When downscaling predictions for several locations in a landscape, spatial correlation of predictions is also important.
With the aim to assist research of climatic impacts on forests, the R package meteoland provides utilities to estimate daily weather variables at any position over complex terrains:
- Spatial interpolation of daily weather records from meteorological stations. }
- Statistical correction of meteorological data series (e.g. from climate models).}
Using meteoland package is easy, but some ideas and concepts must be addresed
to make the experience easy for new users. This vignette is intended as a working example to explain all the steps needed to get comfortable with meteoland workflow. We assume that the reader is more or less familiar with S4 classes of package sp. You will find detailed information regarding the package structure, its functions, and how calculations are done in vignette "Package meteoland".
Installing and loading the package
First of all, before starting to work with the package, we must install and load
the library. meteoland stable version can be found at CRAN
(https://CRAN.R-project.org/package=meteoland), and it can be
installed and loaded as any other R package:
install.packages("meteoland") library(meteoland)
Alternatively, you can install the development version located at GitHub using the devtools package:
library(devtools) install_github("miquelcaceres/meteoland", build_vignettes = TRUE) library(meteoland)
Interpolation of daily weather
Spatial interpolation is required when meteorology for the area and period of interest cannot be obtained from local sensors. The nearest weather station may not have data for the period of interest or it may be located too far away to be representative of the target area.
Preparing weather data for interpolation
data("exampleinterpolationdata") st_names = row.names(exampleinterpolationdata@MinTemperature) st_coords = exampleinterpolationdata@coords st_data = data.frame(X_UTM = st_coords[,1], Y_UTM = st_coords[,2], elevation = exampleinterpolationdata@elevation, row.names = st_names) tmin = exampleinterpolationdata@MinTemperature tmax = exampleinterpolationdata@MaxTemperature prec = exampleinterpolationdata@Precipitation rhum = exampleinterpolationdata@RelativeHumidity
Before starting using the package, you need to have access to the elevation (in m) and daily weather data series corresponding a set of locations (normally weather stations). Elevation is needed because interpolation routines perform corrections for differences in elevation between the reference locations and the target point. The initial format of your data will be different depending on the format used by your data provider (the package has also tools to access weather data). For our example, we will assume you have data from a set of 38 stations in your study area. On one hand, you should have a data.frame with the coordinates and elevation of each location:
str(st_data) head(st_data)
On the other, you should have at least three matrices of meteorological data (one for minimum temperature, one for maximum temperature and the last one for precipitation) with stations in rows and dates in columns. In our example we also add relative humidity (in percent), so that other derived variables can be calculated:
dim(tmax) dim(tmin) dim(prec) dim(rhum) tmax[1:6,1:6]
Units should be in degrees Celsius for temperature and mm for precipitation.
Building an interpolation data object
Meteoland stores weather series for reference locations and interpolation parameters in a single object of class MeteorologyInterpolationData. There are several ways of building such objects, but we will first illustrate how to do it from the data we just presented.
First we need to create an object of class SpatialPoints (see package sp) with the spatial coordinates of our stations and the coordinate system (here UTM 31N):
sp = SpatialPoints(st_data[,c("X_UTM", "Y_UTM")], proj4string = CRS("+proj=utm +zone=31 +ellps=WGS84 +datum=WGS84 +units=m +towgs84=0,0,0")) head(sp)
We can now build an object MeteorologyInterpolationData using:
interpolator <- MeteorologyInterpolationData(sp, elevation = st_data$elevation, MinTemperature = tmin, MaxTemperature = tmax, Precipitation = prec, RelativeHumidity = rhum) class(interpolator)
The resulting object is now ready to be used to perform interpolation on a set of target locations (see next section). We can inspect the amount of data in our interpolation object using function interpolation.coverage. For example, with:
spatial_coverage <- interpolation.coverage(interpolator, type = 'spatial') head(spatial_coverage)
we obtain the number of non-missing observations for each weather station and variable. Similarly, we can ask the number of non-missing observations for each date and variable using:
temporal_coverage <- interpolation.coverage(interpolator, type = 'temporal') head(temporal_coverage)
Interpolation parameters are also stored in the same object (see next subsection):
names(interpolator@params)
Interpolation basics and parameters
Package meteoland implements, with a few modifications, the daily weather interpolation and estimation algorithms that un- derpin the U.S. DAYMET dataset (Thornton et al., 1997; Thornton and Running, 1999).This approach, similar to inverse distance weighting, interpolates weather variables using trun- cated Gaussian filters, which consist in defining spatial weights $W(r)$ at radial distance $r$ from a target point $p$ using: \begin{equation} W(r) = e^{-\alpha \cdot (r/R_p)^2} - e^{-\alpha} \end{equation} if $r \leq R_p$ and $W(r) = 0$ otherwise.
Here $r$ is the radial distance from $p$, $R_p$ is the truncation distance and $\alpha$ is the shape parameter. The spatial convolution of this filter with a set of weather station locations results, for each target point, in a vector of weights associated with observations. The following figure illustrates the Gaussian filter for $R_p = 500$ and either $\alpha = 3.0$ (continuous line) or $\alpha = 6.25$ (dashed line):
r = 0:1000 R_p = 500 gf1 = exp(-3.0*((r/R_p)^2.0)) - exp(-3.0) gf2 = exp(-6.25*((r/R_p)^2.0)) - exp(-6.25) gf1[r>R_p] = 0 gf2[r>R_p] = 0 plot(r, gf1, type="l", ylab = "W(r)", xlab ="r") lines(r, gf2, lty=2) legend("topright", legend = c("alpha = 3", "alpha = 6.25"), lty=c(1,2), bty="n")
$R_p$ is automatically adjusted so that it has lower values in data-rich regions and is increased in data-poor regions. The method, however, requires the user to specify $N$, the average number of observations to be included for each target point. $R_p$ is then varied as a smooth function of the local density in such a way that this average is achieved over the spatial domain. It is important that the initial value of $R_p$
In meteoland estimation of $Rp$ is done once for each target point, variable and day. Interpolation of temperature includes a correction for the effects of elevation. More specifically, a weighted least-squares regression is used to assess the relationship between temperature differences and elevation differences in weather station data and this relationship is applied to elevation differences between weather stations and the target point. Interpolation of relative humidity is done after transforming it to dew-point temperature. No correction for elevation is performed during interpolation, but elevation effects arise when back-transforming dew-point temperature to relative humidity. Interpolation of daily precipitation is complicated by the need to predict both precipitation occurrence and, conditioned on this, precipitation amount. Thornton et al. (1997) defined a binomial predictor of spatial precipitation occurrence as a function of the weighted occurrence at surrounding weather stations. Conditional on precipitation occurrence, the interpolation routine predicts precipitation amount, where weighted least-squares regression is also used to account for elevation effects. Interpolation of wind is performed in three different ways depending on the information available. If only wind speed data is available, the spatial interpolation with Gaussian weights is used on wind scalars as described above. If weather station data includes wind direction, a polar average is cal- culated using Gaussian weights. Finally, if static wind fields are also available the interpolation routine first finds, for each weather station the wind field that best matches the observed vector. Then, the wind vectors extracted from the selected wind fields are averaged as before. The details of how $R_p$ is estimated and, in general, how interpolation is done are given in a PDF vignette called 'Meteorology'.
Interpolation parameters $\alpha$ and $N$ can be different for each variable to be interpolated. The following table lists all the interpolation parameters (see also function defaultInterpolationParameters())
| Paremeter | Default value | Definition |
| ------------- | ---------- | ----------------------------------- |
| initial_Rp | 140000 | Initial truncation radius |
| iterations | 3 | Number of station density iterations |
| alpha_MinTemperature | 3.0 | Gaussian shape parameter for minimum temperature |
| alpha_MaxTemperature | 3.0 | Gaussian shape parameter for maximum temperature |
| alpha_DewTemperature | 3.0 | Gaussian shape parameter for dew-point temperature |
| alpha_PrecipitationEvent | 5.0 | Gaussian shape parameter for precipitation events |
| alpha_PrecipitationAmount | 5.0 | Gaussian shape parameter for the regression of precipitation amounts |
| alpha_Wind | 3.0 | Gaussian shape parameter for wind |
| N_MinTemperature | 30 | Average number of stations with non-zero weights for minimum temperature |
| N_MaxTemperature | 30 | Average number of stations with non-zero weights for maximum temperature |
| N_DewTemperature | 30 | Average number of stations with non-zero weights for dew-point temperature |
| N_PrecipitationEvent | 5 | Average number of stations with non-zero weights for precipitation events |
| N_PrecipitationAmount | 20 | Average number of stations with non-zero weights for the regression of precipitation amounts |
| N_Wind | 2 | Average number of stations with non-zero weights for wind |
| St_Precipitation | 5 | Number of days for the temporal smoothing of precipitation |
| St_TemperatureRange | 15 | Number of days for the temporal smoothing of temperature range |
| pop_crit | 0.50 | Critical precipitation occurrence parameter |
| f_max | 0.6 | Maximum value for precipitation regression extrapolations (0.6 equals to a maximum of 4 times extrapolation) |
| wind_height | 10 | Wind measurement height (in m) |
Parameter St_Precipitation controls the temporal smoothing that is applied to weather station data to calibrate regression, while parameters pop_crit and f_max are also particular to the estimation of precipitation (see details in Thornton et al. 1997). Parameter St_TemperatureRange is used for the estimation of solar radiation (see details in PDF vignette called 'Meteorology').
A parameter that is particularly important to understand is initial_Rp, which specifies the initial radius for the truncated spatial Gaussian kernel. By default its value is:
interpolator@params$initial_Rp
The value of initial_Rp must be set in relation to the units of the spatial coordinates of weather data. Our data was in meters, so the default radius is 140 km. In general, the initial radius should be large enough to include a reasonable number of stations (~20-40), but the kernel radius is adjusted for each interpolation target point.
Calibration and cross-validation of the interpolation data
Once we already have weather stations data in shape, we can start calibrating the model in order to obtain the optimal parameters for the meteorological variables we want to interpolate. Parameter calibration has to be done for each variable separately. For example for minimum temperature:
tmin_cal <- interpolation.calibration(interpolator, variable = "Tmin", N_seq = 20, alpha_seq = seq(5, 10, by = 1), verbose = TRUE)
This function returns an interpolation.calibration class object which contains several items:
- Numeric matrix with the mean absolute error (MAE) values for each combination of parameters $N$ and $\alpha$.
- Miminum value found for MAE.
- Value for the $N$ parameter corresponding to the minumun MAE.
- Value for the $\alpha$ parameter corresponding to the minimum MAE.
- Matrix with the observed values.
- Matrix with the predicted values for the optimum parameter combination.
The result of the calibration needs to be manually stored in the interpolation params:
interpolator@params$N_MinTemperature = tmin_cal$N interpolator@params$alpha_MinTemperature = tmin_cal$alpha
We strongly recommend conducting calibration exercises at least once for each variable and each data set used as reference for interpolation, and more than once if periods differ in the number of stations available.
Before using the object for interpolations, we also need to assess its performance. This is done by cross-validation in function interpolation.cv:
cv <- interpolation.cv(interpolator, verbose=T)
Cross-validation is perfomed by leave-one-out, which means that each target station is first removed from the data set and then the remaining stations are used to obtain interpolation estimates. In this way one maximizes the amount of information for estimates, while keeping them independent of the observed values in each station. The results of cross-validation can be inspected using a specific summary function:
summary(cv)
Interpolation on a grid
data("examplegridtopography") elev = as(examplegridtopography[1:20,1:20], "SpatialGridDataFrame")["elevation"]
The target for weather interpolation in meteoland can be a set of points, pixels or a whole grid. Again, the initial format of data can be very different. Here we assume you have a small grid of 400 (20x20) cells of 1ha in size, with elevation data in form of class SpatialGridDataFrame (see method read.asciigrid in package sp):
summary(elev)
Note that the coordinate reference system needs to be the same as that of interpolator, which in this case it is. Before performing the interpolation over this grid, we need to reshape this data in a class called SpatialGridTopography:
sgt = SpatialGridTopography(as(elev, "SpatialGrid"), elevation = elev$elevation, proj4string = elev@proj4string) sgt
As you can see in the result, meteoland has calculated for us slope and aspect (both in degrees) from elevation data ^[Objects of class SpatialGridTopography can be initialized with user input values for slope and aspect too, but meteoland has its own routines when this are missing]. Slope and aspect are important for radiation calculations, which also requires relative humidity data. We can display elevation over the grid using:
spplot(sgt, "elevation")
Before we call the interpolation routine, we need to define the dates (i.e. days) for which we want weather to be interpolated, for example:
dates = as.Date(c("2001-02-03", "2001-06-03"))
Of course, we need to be sure that the interpolator object has data corresponding to this dates. We can check if there is any missing date using:
sum(!(dates %in% interpolator@dates))
The name of interpolation functions depend on the target spatial structure. For grids we need to use function interpolationgrid:
ml <- interpolationgrid(interpolator, sgt, dates)
This function works processing each date at a time. Since calculations can take some time, the console output shows the progress. The output of the function is an object of class SpatialGridMeteorology:
ml
We can display interpolated grids in a map using function spplot:
spplot(ml, 2, "MinTemperature") spplot(ml, 2, "MaxTemperature")
Objects of class SpatialGridMeteorology include a list of data frames, one per date. We can access the interpolated data for a given date using:
df_1 = ml@data[[1]] head(df_1)
Some columns are missing (e.g. wind speed) because we did not include weather station data regarding these variables. If we wants to add grid coordinates to this data frame, we can use the spatial information stored in the ml object:
sgdf_1 = SpatialGridDataFrame(grid = ml@grid, data = ml@data[[1]], proj4string = ml@proj4string) summary(sgdf_1)
Results can be retrieved in this way and saved using write.asciigrid for its use outside R. The package also provides function for reading/writing NetCDFs (see function writemeteorologygrid) from SpatialGridMeteorology objects.
Interpolation on a set of points
spt = as(examplegridtopography, "SpatialPointsTopography")[c(36,228, 530, 2304)] cc = coordinates(spt) points_df = data.frame(X_UTM = cc[,1], Y_UTM = cc[,2], elevation = spt@data$elevation)
If you want to interpolate on a set of target locations, the starting point will normally be a data.frame. In our example this points come from the grid, but we have reshaped them so the starting format is familiar:
points_df
Analogously with the grid, we need to transform this data into an object SpatialPointsTopography:
spt = SpatialPointsTopography(as.matrix(points_df[,c("X_UTM", "Y_UTM")]), elevation = points_df$elevation, proj4string = CRS("+proj=utm +zone=31 +ellps=WGS84 +datum=WGS84 +units=m +towgs84=0,0,0")) spt
In this case we only have elevation (in m), but slope and aspect should also be included if possible. Let us assume you want to interpolate on this points for the whole time series available in object interpolator. Since we are dealing with points, the function to interpolate is called interpolationpoints:
mp = interpolationpoints(interpolator, spt)
This function works processing one point at a time. The output of the function is an object of class SpatialPointsMeteorology:
mp
And the time series for a given can be plotted using function meteoplot. For example, we show here the precipitation series of point #1:
meteoplot(mp, 1, "Precipitation", ylab="Precipitation (mm)", xlab="")
Objects of class SpatialPointsDataFrame include a list of data frames, one per point. We can access one of them using:
df_1 = mp@data[[1]] head(df_1)
This data frame can now be written into a file for its analysis outside R. The package also provides its own functions to write/read point meteorology data in different formats. If we are interested in inspecting the interpolation result by date, instead of by point, we can use function extractpointdates, which returns objects of class SpatialGridDataFrame:
dt_4 = extractpointdates(mp, as.Date("2001-01-04"), verbose = FALSE) dt_4
Statistical correction of daily weather
Correcting the biases of a meteorological data series containing biases using a more accurate meteorological series is necessary when the more accurate series does not cover the period of interest and the less accurate series does. The less accurate series may be at coarser scale, as with climate model predictions or climate reanalysis data. In this case one can speak of statistical correction and downscaling. However, one may also correct the predictions of climate models using reanalysis data estimated at the same spatial resolution.
In the following example we will correct the predictions of Regional Climate Model (CCLM4-8-17; driving global model CNRM-CERFACS-CNRM-CM5) on the same area of the interpolation example. RCM data includes 3 model cells. Meteorological data covers an historical (reference) period (2000-2003) and a future (projection) period (year 2023), the latter simulated under rcp4.5 scenario.
Preparing weather data to be downscaled/corrected
data("examplecorrectiondata") coords = examplecorrectiondata@coords pt_coords = data.frame(long = as.numeric(coords[,1]), lat = as.numeric(coords[,2])) proj_dates = examplecorrectiondata@dates ref_data = examplecorrectiondata@reference_data proj_data = examplecorrectiondata@projection_data
One needs several data items to perform downscaling and statistical correction. First, we need a matrix or data frame with the central coordinates of the RCM cells (here in longitude/latitude format):
pt_coords
Second, we need the uncorrected meteorological data (here RCM outputs) for both a reference period (that will be matched with our more accurate series) and a projection period (that our more accurate series does not cover). Both should be arranged in lists of data frames, one per RCM cell:
length(ref_data) length(proj_data) head(ref_data[[1]]) head(proj_data[[1]])
Note that these data frames should follow the conventions of meteoland for variable names and units. Finally, we need to specify the projection period, in this case one year:
proj_dates = seq(as.Date("2023-01-01"), as.Date("2023-12-31"), by="day")
Building the uncorrected data object
Statistical correction needs an object of class MeteorologyUncorrectedData, analogous to the MeteorologyInterpolationData of the previous section. To build this object we need first to express coordinates in an object SpatialPoints:
sp = SpatialPoints(pt_coords, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")) sp
Assuming the weather data to be corrected is in the proper format, the object is constructed using:
uncorrected = MeteorologyUncorrectedData(sp, ref_data, proj_data, proj_dates)
Target point weather meteorology taken as reference
Downscaling/statistical correction is done on a set of target spatial points. meteoland will first find the RCM cell to which each point is nearest (using the central coordinates of RCM cells). Once this matching is done, statistical relations are build between the meteorology of RCM cells and that of target points for the reference period (here 2000-2003) which are used to correct RCM data for the projection period (year 2023). For this process to be done, we need the coordinates and weather series for the reference period. In our case, we will employ a copy the object of class SpatialPointsMeteorology that resulted from interpolation in the previous section:
historical = mp
historical
We named it historical to remind that we are using the historical period as reference and want to correct the RCM series for the projection period There are several ways to shape data into objects of class SpatialPointsMeteorology. For example, one can use function readmeteorologypoints.
Tuning correction parameters
Objects of class MeteorologyUncorrectedData have a slot with parameters:
uncorrected@params
Importantly, varmethods is a named list that specifies which correction method has to be used for each variable. Since we do not have reference data for wind speed, we must turn the method for wind speed to "none":
uncorrected@params$varmethods$WindSpeed="none"
Conducting statistical correction
Once we have all the data objects that we need (the hard part), conducting statistical correction is straightforward:
projected = correctionpoints(uncorrected, historical)
As mentioned above, correction proceeds point by point. The result is again an object of class SpatialPointsMeteorology:
projected
The following code displays the minimum/maximum temperatures before and after correction:
#Plot predicted mean temperature for point 1 meteoplot(projected, 1, "MinTemperature", ylab="Temperature (Celsius)", ylim=c(-5,40), col="blue") meteoplot(projected, 1, "MaxTemperature", add=TRUE, col="red") #Add uncorrected mean temperature data (cell #3) lines(uncorrected@dates, uncorrected@projection_data[[3]]$MinTemperature, col="blue", lty=3) lines(uncorrected@dates, uncorrected@projection_data[[3]]$MaxTemperature, col="red", lty=3) legend("topright", legend=c("corrected","uncorrected", "Maximum", "Minimum"), col=c("black","black", "red","blue"), lty=c(1,3,1,1), bty="n")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.