Skip to content

Regridding and temporal aggregation (data collocation)

Jump to bottom
Sixto Herrera García edited this page Feb 19, 2016 · 4 revisions

Regridding and interpolation

Function interpData in loadeR performs interpolation of a gridded or points dataset into a new user-defined spatial pattern using bilinear weights or nearest-neighbour methods. In the following, we illustrate the different options that are supported using gridded data from the NCEP/NCAR Reanalysis 1 dataset and point (station) data from the VALUE ECA&D observational dataset (see previous section for more information about these datasets).

Grid data:

# Download NCEP (model data) datasets
dir.create("mydirectory")
download.file("http://meteo.unican.es/work/loadeR/data/Iberia_NCEP.tar.gz",
              destfile = "mydirectory/Iberia_NCEP.tar.gz")
# Extract files from the tar.gz file
untar("mydirectory/NCEP_Iberia.tar.gz", exdir = "mydirectory")
# Path to the NCEP ncml file.
ncep <- "mydirectory/Iberia_NCEP/Iberia_NCEP.ncml"

# Load air temperature at 1000 mb isobaric pressure level for boreal winter (DJF) 1991-2000
t1000.djf <- loadGridData(ncep, 
                          var = "ta@1000", 
                          lonLim = c(-12,10), 
                          latLim = c(33,47), 
                          season = c(12,1,2), 
                          years = 1991:2000,
                          dictionary = TRUE)

## [2016-02-18 10:47:27] Defining homogeneization parameters for variable "ta"
## [2016-02-18 10:47:27] Defining geo-location parameters
## [2016-02-18 10:47:27] Defining time selection parameters
## [2016-02-18 10:47:28] Retrieving data subset ...
## [2016-02-18 10:47:28] Done
plotMeanField(t1000.djf)

Station data:

download.file("http://meteo.unican.es/work/loadeR/data/VALUE_ECA_86_v2.tar.gz", 
              destfile = "mydirectory/VALUE_ECA_86_v2.tar.gz")
# Extract files from the tar.gz file
untar("mydirectory/VALUE_ECA_86_v2.tar.gz", exdir = "mydirectory")
# Data inventory
value <- "mydirectory/VALUE_ECA_86_v2"

obs <- loadStationData(dataset = value, var= "tmean", lonLim = c(-12,10), latLim = c(33,47), 
                       season = c(12,1,2), years = 1991:2000)

## [2016-02-18 10:47:30] Loading data ...
## [2016-02-18 10:47:30] Retrieving metadata ...
## [2016-02-18 10:47:30] Done.

x <- getCoordinates(obs)$x
y <- getCoordinates(obs)$y
plot(x,y, pch = 10, col = "blue")

library(fields)
world(add=TRUE)

From grid to grid

# Bilinear interpolation to a smaller domain centered in Spain using a 0.5 degree resolution
# in both X and Y axes
t1000.djf.05 <- interpData(t1000.djf, 
                           new.Coordinates = list(x = c(-10,5,.5), y = c(36,44,.5)),
                           method = "bilinear")

## [2016-02-18 10:47:31] Performing bilinear interpolation... may take a while
## [2016-02-18 10:47:31] Done

plotMeanField(t1000.djf.05)

Note that new attributes "interpolation", "resX" and "resY" indicate that the original data have been interpolated:

attributes(t1000.djf.05$xyCoords)

## $names
## [1] "x" "y"
## 
## $projection
## [1] "LatLonProjection"
## 
## $resX
## [1] 0.5
## 
## $resY
## [1] 0.5
## 
## $interpolation
## [1] "bilinear"

We could use method = "nearest" to perform the nearest neighbour interpolation in all the cases of this worked example, for instance:

t1000.djf.05 <- interpData(t1000.djf, 
                           new.Coordinates = list(x = c(-10,5,.5), y = c(36,44,.5)),
                           method = "nearest")



## [2016-02-18 10:47:32] Performing nearest interpolation... may take a while
## [2016-02-18 10:47:44] Done

plotMeanField(t1000.djf.05)

The new.Coordinatesargument can be filled with an existing grid of other loaded data. We can use function getGridto extract the coordinates:

data(tp_forecast)
par(mfrow = c(1,3))
plotMeanField(tp_forecast)
plotMeanField(t1000.djf)

# Bilinear interpolation to a smaller domain centered in Spain using a 0.5 degree resolution
# in both X and Y axes
new.Coordinates <- getGrid(tp_forecast)
t1000.djf.new <- interpData(t1000.djf, new.Coordinates = new.Coordinates,
 method = "bilinear")

## [2016-02-18 10:47:45] Performing bilinear interpolation... may take a while
## [2016-02-18 10:47:45] Done

plotMeanField(t1000.djf.new)

From grid to points

new.Coordinates <- list(x = obs$xyCoords[,1], y = obs$xyCoords[,2])
attr(new.Coordinates,"type") <- "location"

t1000.djf.new <- interpData(t1000.djf, new.Coordinates = new.Coordinates, method = "bilinear")

## [2016-02-18 10:47:46] Performing bilinear interpolation... may take a while
## [2016-02-18 10:47:47] Done

x <- getCoordinates(t1000.djf.new)$x
y <- getCoordinates(t1000.djf.new)$y
plot(x,y,  pch = 10, col = "blue")

library(fields)
world(add=TRUE)

From points to grid

new.Coordinates <- list(x = t1000.djf$xyCoords$x, y = t1000.djf$xyCoords$y)
# or...
new.Coordinates <- getGrid(t1000.djf)

t1000.djf.new <- interpData(obs, new.Coordinates = new.Coordinates, method = "bilinear")

plotMeanField(t1000.djf.new)

Time aggregation

Function timeAggregation in loadeR performs data transformation to the desired time scale ("daily", "monthly" or "annual") by specifying an aggregation function. In this example, daily precipitation data is loaded and the annual accumulated precipitation (sum of the daily precipitation) is obtained:

  • Download NCEP (model data) datasets:
dir.create("mydirectory")
download.file("http://meteo.unican.es/work/loadeR/data/Iberia_NCEP.tar.gz",
              destfile = "mydirectory/Iberia_NCEP.tar.gz")
# Extract files from the tar.gz file
untar("mydirectory/NCEP_Iberia.tar.gz", exdir = "mydirectory")
# Path to the NCEP ncml file.
ncep <- "mydirectory/Iberia_NCEP/Iberia_NCEP.ncml"
  • Loading precipitation for the period 1991-2000:
tp <- loadGridData(ncep, 
                   var = "tp", 
                   lonLim = c(-12,10), 
                   latLim = c(33,47),
                   years = 1991:2000,
                   dictionary = TRUE)

## [2016-02-18 10:47:47] Defining homogeneization parameters for variable "tp"
## [2016-02-18 10:47:47] Defining geo-location parameters
## [2016-02-18 10:47:47] Defining time selection parameters
## [2016-02-18 10:47:48] Retrieving data subset ...
## [2016-02-18 10:47:48] Done
  • Obtaining the annual accumulated precipitation:
tp_annual <- timeAggregation(tp, aggr.d = "sum", aggr.m = "sum", aggr.y = "sum")

## Data is not sub-daily: aggr.d will be ignored
## Applying sum monthly aggregation function
## [2016-02-18 10:47:48] Aggregating...
## [2016-02-18 10:47:48] Done.
## Applying sum annual aggregation function
## [2016-02-18 10:47:48] Aggregating...
## [2016-02-18 10:47:48] Done.

The message returned by the function reports that the input data is not sub-dalily, thus we could omit the aggr.d argument:

tp_annual <- timeAggregation(tp, aggr.m = "sum", aggr.y = "sum")

## Applying sum monthly aggregation function
## [2016-02-18 10:47:48] Aggregating...
## [2016-02-18 10:47:48] Done.
## Applying sum annual aggregation function
## [2016-02-18 10:47:48] Aggregating...
## [2016-02-18 10:47:48] Done.

Note that the aggregation is done by steps, this is, data is first monthly aggregated and then annually. However, timeAggregation allows the direct aggregation of daily data into annual data:

tp_annual <- timeAggregation(tp, aggr.y = "sum")

## Applying sum annual aggregation function
## [2016-02-18 10:47:49] Aggregating...
## [2016-02-18 10:47:49] Done.

To obtain monthly data, argument aggr.y is omitted:

tp_monthly <- timeAggregation(tp, aggr.m = "sum")

## Applying sum monthly aggregation function
## [2016-02-18 10:47:49] Aggregating...
## [2016-02-18 10:47:49] Done.

Monthly data can be also the input:

tp_annual <- timeAggregation(tp_monthly, aggr.y = "sum")

## Applying sum annual aggregation function
## [2016-02-18 10:47:49] Aggregating...
## [2016-02-18 10:47:49] Done.
# Plots of the daily and annual precipitation
par(mfrow = c(1,2))
plotMeanField(tp)
plotMeanField(tp_annual)

<-- Home page of the Wiki

print(sessionInfo())

## R version 3.2.3 (2015-12-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.3 LTS
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
##  [3] LC_TIME=es_ES.UTF-8           LC_COLLATE=en_US.UTF-8       
##  [5] LC_MONETARY=es_ES.UTF-8       LC_MESSAGES=en_US.UTF-8      
##  [7] LC_PAPER=es_ES.UTF-8          LC_NAME=es_ES.UTF-8          
##  [9] LC_ADDRESS=es_ES.UTF-8        LC_TELEPHONE=es_ES.UTF-8     
## [11] LC_MEASUREMENT=es_ES.UTF-8    LC_IDENTIFICATION=es_ES.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] fields_8.3-6      maps_3.0.2        spam_1.3-0        loadeR_0.2-0     
## [5] loadeR.java_1.1-0 rJava_0.9-8      
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-31 digest_0.6.8    bitops_1.0-6    formatR_1.2    
##  [5] magrittr_1.5    evaluate_0.7    stringi_0.4-1   sp_1.1-0       
##  [9] akima_0.5-12    rmarkdown_0.6.1 tools_3.2.3     stringr_1.0.0  
## [13] RCurl_1.95-4.7  yaml_2.1.13     abind_1.4-3     htmltools_0.2.6
## [17] knitr_1.10.5

loadeR - Santander MetGroup (Univ. Cantabria - CSIC)

  1. Package Installation (and known problems)
  2. Model Data (reanalysis and climate projections)
    1. Dataset definition and loading local grid data
    2. Loading remote grid data
    3. Loading data from the User Data Gateway (UDG)
  3. Observations (station and gridded data)
    1. Standard (ASCII) format for station data
    2. Accessing and loading station data
    3. Accessing and loading grid data
  4. Standard data manipulation
    1. Harmonization
    2. Exporting to netCDF
Clone this wiki locally