Changes between Version 41 and Version 42 of udg/ecoms/RPackage/examples


Ignore:
Timestamp:
Feb 14, 2014 3:00:37 PM (8 years ago)
Author:
juaco
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • udg/ecoms/RPackage/examples

    v41 v42  
     1# EXAMPLE: Loading and plotting various members:
    12
    2 The following examples are partly based on locally stored data. Only the `loadSeasonalForecast` function does not require locally stored datasets because it works by remotely accessing the System4 databases stored in the SPECS-EUPORIAS Data Server (in this case a valid username and password are required though, as detailed [https://www.meteo.unican.es/trac/meteo/wiki/EcomsUdg/RPackage/Authentication Authentication section]). Therefore, it is recommended that the [https://www.meteo.unican.es/trac/meteo/raw-attachment/wiki/EcomsUdg/meteoR_v1_0.zip meteoR_v1_0.zip] file is downloaded and unzipped in a convenient directory before starting. Please follow the [wiki:../Requirements instructions].
     3# Total precipitation of January 2010 forecasted in October 2009 (lead month 3) by the System4 model (seasonal range, 15 members) 
     4# is next represented for each member, using the spplot method for the SPatialGridDataFrame class of library sp:
     5
     6gg.pr <- loadSeasonalForecast("System4_seasonal_15", var="tp", members=NULL, lonLim=c(-30,20), latLim=c(-12,15), season=1, years=2010, leadMonth=3)
     7
     8names(gg.pr)
     9str(gg.pr$MemberData)
    310
    411
    5 = loadSeasonalForecast
     12pr.list <- lapply(1:length(gg.pr$MemberData), function(x) colSums(gg.pr$MemberData[[x]]))
     13df <- do.call("data.frame", pr.list)
     14names(df) <- names(gg.pr$MemberData)
    615
    7 In the next lines we describe an illustrative example of the `loadSeasonalForecast` function. '''Note that a valid username and password are required''' in order to access the SPECS-EUPORIAS Data Server, as detailed [https://www.meteo.unican.es/trac/meteo/wiki/EcomsUdg/RPackage/Authentication here]. We will retrieve System4 simulation data for the Iberian Peninsula, considering mean surface temperature for January and the first simulation member, for the 10-year period 1990-1999. This simple example has been chosen because of the fast data access (note that this also depends on the connection speed). Using a standard broadband connection, loading the following dataset took approximately 19 seconds:
     16class(gg.pr$LonLatCoords)
     17sgdf <- SpatialGridDataFrame(gg.pr$LonLatCoords, df)
    818
    9 {{{
    10 > openDAP.query <- loadSeasonalForecast(dataset = "http://www.meteo.unican.es/tds5/dodsC/system4/System4_Seasonal_15Members.ncml",
    11 +                                         standard.vars = TRUE, dictionary = "datasets/forecasts/System4/System4_Seasonal_15Members.dic",
    12 +                                         var = "tas", members = 1,
    13 +                                         lonLim = c(-10,5), latLim = c(35,45),
    14 +                                         season = 1, years = 1990:1999, leadMonth = 1)
    15 
    16 }}}
     19spplot(sgdf, scales=list(draw = TRUE), col.regions = rev(terrain.colors(50)), at = seq(0, ceiling(max(sgdf@data)),10))
    1720
    1821
    19 Data are now ready for analysis into our R session. Note that we have set the `standard.vars` argument to `TRUE`, and we have specified a path to the dictionary, a file with extension ''.dic'', whose aim is to translate the original variable stored in the dataset (in this case in Kelvin instead of Celsius). More details on the use of the dictionary are provided in [wiki:EcomsUdg/RPackage this section].
     22data(world_map)
     23wl <- as(world_map, "SpatialLines")
     24l1 <- list("sp.lines", wl)
    2025
    21 {{{
    22 > str(openDAP.query)
    23 List of 6
    24  $ VarName      : chr "tas"
    25  $ isStandard   : logi TRUE
    26  $ MemberData   :List of 1
    27   ..$ : num [1:310, 1:280] 13.3 13.9 12.5 13 13 ...
    28  $ LatLonCoords : num [1:280, 1:2] 45 44.2 43.5 42.7 42 ...
    29   ..- attr(*, "dimnames")=List of 2
    30   .. ..$ : NULL
    31   .. ..$ : chr [1:2] "lat" "lon"
    32  $ RunDates     : chr [1:310] "1989-12-01" "1989-12-01" "1989-12-01" "1989-12-01" ...
    33  $ ForecastDates:List of 2
    34   ..$ Start: POSIXlt[1:310], format: "1990-01-01" "1990-01-02" "1990-01-03" "1990-01-04" ...
    35   ..$ End  : POSIXlt[1:310], format: "1990-01-02" "1990-01-03" "1990-01-04" "1990-01-05" ...
    36 
    37 }}}
    38 
    39 
    40 The function has been implemented to access seasonal slices (as determined by the `season` argument. Seasons can be defined in several ways: A single month (e.g. `season = 1` for January, as in this example), a standard season (e.g. `season=1:3` for JFM, or `season=c(12,1,2)` for DJF), or any period of consecutive months (e.g. `season=1:6`, for the first half of the year). Seasons are returned for a given year period (defined by the `years` argument, e.g. `years = 1990:1999` in this example) with a homogeneous forecast lead time (as given by the `leadMonth` argument; in this example `leadMonth = 1` for one-month lead time) with respect to the first month of the selected season. For example, in this particular case the data loaded correspond to the series of January 1990 to January 1999 from the December 1989 to December 1998 runtime forecast. As a result, the length of the time series returned for each of the 280 grid cells is of 310 days (31 days of January * 10 years). Note that it is also possible to work with year-crossing seasons, such as DJF (in this case `season=c(12,1,2)`).
    41 
    42 
    43 
    44 A common task consists of the representation of data, e.g. by mapping the spatial mean for the period considered. Another common task is the representation of time series for selected point locations/grid cells. In this example, we will map the mean temperature field for the period selected (1990-99) preserving the original spatial resolution of the model. Furthermore, we will display time series at two grid points coincident with the locations of two Spanish cities. To this aim, we will make use of some `base` R functions and also from some contributed packages that can be very useful for climate data handling and representation.
    45 
    46 {{{
    47 > # Calculation of the mean values of the period for each grid cell:
    48 > mean.field <- colMeans(openDAP.query$MemberData[[1]])
    49 > # Creation of a matrix with selected point locations:
    50 > city.names <- c("Madrid", "Santander")
    51 > locations <- matrix(c(-3.68, 40.4, -3.817, 43.43), ncol=2, byrow = TRUE)
    52 > dimnames(locations) <- list(city.names, c("lon","lat"))
    53 > print(locations)
    54              lon   lat
    55 Madrid    -3.680 40.40
    56 Santander -3.817 43.43
    57 
    58 }}}
    59 
    60 Note that the geographical coordinates of the requested spatial domain are nor perfectly uniform, and as a result it is not possible to represent the data as a regular grid. To overcome this problem, we perform an interpolation, although we preserve the native grid cell size of the model (0.75deg) for data representation. The R package [http://cran.r-project.org/web/packages/akima/index.html akima] provides a extremely fast interpolation algorithm by means of the `interp` function.
    61 
    62 {{{
    63 > lat <- openDAP.query$LatLonCoords[ ,1]
    64 > lon <- openDAP.query$LatLonCoords[ ,2]
    65 # Requires "akima::interp" for regular grid (bilinear) interpolation
    66 > library(akima)
    67 > grid.075 <- interp(lon, lat, mean.field, xo = seq(min(lon), max(lon), .75), yo = seq(min(lat), max(lat), .75))
    68 
    69 }}}
    70 
    71 After the interpolation, data are now in a regular grid. Note that the object `grid.075` is a list with the usual `x,y,z` components as required by many R functions for gridded data representation (e.g. `image`, `contour`, `persp`...)
    72 
    73 {{{
    74 > str(grid.075)
    75 List of 3
    76  $ x: num [1:19] -9.75 -9 -8.25 -7.5 -6.75 ...
    77  $ y: num [1:14] 35.2 36 36.7 37.5 38.2 ...
    78  $ z: num [1:19, 1:14] 15.1 15.1 15 15 14.3 ...
    79 
    80 }}}
    81 
    82 In the following lines of code we plot the mean temperature field. In addition, we will also add to the map the point locations of the selected cities for which the time series will be represented. The R package [http://cran.r-project.org/web/packages/fields/index.html fields] provides many useful tools for spatial data handling and representation, including a world map that can be easily incorporated in our plots.
    83 
    84 {{{
    85 > # Representation of the mean temperature of the period
    86 > library(fields)
    87 > image.plot(grid.075, asp=1,ylab = "latitude", xlab = "longitude", col=topo.colors(36),
    88 +           main = "Mean surface T January 1990-99", legend.args=list(text="degC", side=3, line=1))
    89 # Adds selected locations to the plot and puts labels:
    90 > points(locations, pch=15)
    91 > text(locations, pos=3, city.names)
    92 # Adds the world map to the current plot:
    93 > world(add=TRUE)
    94 
    95 }}}
    96 
    97 [[Image(TmeanJanS4.png)]]
    98 
    99 Next, we plot the time series for the selected locations. To this aim, we calculate the nearest grid points to the specified locations. This can be easily done using the function `fields::rdist`. Note that the output of `loadSeasonalForecast` returns a matrix of Lat-Lon coordinates, as usually found in many climate datasets.  However, the usual format of 2D coordinates matrix in R is Lon-Lat. As a result, note that we specify the coordinates by reversing the column order (i.e.: `openDAP.query$LatLonCoords[ ,2:1]` instead of `openDAP.query$LatLonCoords`):
    100 
    101 {{{
    102 > # Creation of a Euclidean distance matrix among all pairs of selected locations and grid points
    103 > dist.matrix <- rdist(openDAP.query$LatLonCoords[ ,2:1], locations)
    104 > # Positions of the nearest grid points
    105 > index <- rep(NA,ncol(dist.matrix))
    106 > for (i in 1:ncol(dist.matrix)) {
    107 +      index[i] <- which.min(dist.matrix[ ,i])
    108 + }
    109 > # index contains the column positions in the `MemberData` matrices
    110 > locations.data <- openDAP.query$MemberData[[1]][ ,index]
    111 > colnames(locations.data) <- city.names
    112 > str(locations.data)
    113  num [1:310, 1:2] 2.52 3.57 4.05 6.88 7.36 ...
    114  - attr(*, "dimnames")=List of 2
    115   ..$ : NULL
    116   ..$ : chr [1:2] "Madrid" "Santander"
    117 
    118 }}}
    119 
    120 The object `locations.data` is a matrix in which time series are arranged in columns for each of the four locations selected.
    121 
    122 {{{
    123 > ylimits <- c(floor(min(locations.data)), ceiling(max(locations.data)))
    124 > plot(locations.data[ ,1], ty='n', ylim = ylimits, axes=FALSE, ylab="degC", xlab="Year")
    125 > axis(1,at = seq(1,31*11,31), labels=c(1990:1999,""))
    126 > axis(2, ylim=ylimits)
    127 > abline(v=seq(1,31*10,31), lty=2)
    128 > for (i in 1:ncol(locations.data)) {
    129 +      lines(locations.data[ ,i], col=i)
    130 + }
    131 > legend("bottomleft", city.names, lty=1, col=1:4)
    132 > title(main = "Mean surface Temperature January")
    133 > mtext("System4 15 member Seasonal - 1st Member, lead month = 1")
    134 
    135 }}}
    136 
    137 
    138 [[Image(timeSeriesS4.png)]]
    139 
    140 
    141 Alternatively, for selected point locations the function allows for the retrieval of single-point data. In this case, we can enter the lon and lat coordinates of the desired point in the `lonLim`and `latLim` arguments of the function. The function operates by finding the nearest grid point of the dataset to the given coordinates (in terms of Euclidean distances). For instance, the next two instructions load the data represented in the time series above directly for Santander and Madrid respectively:
    142 
    143 {{{
    144 > santanderData <- loadSeasonalForecast(dataset = "http://www.meteo.unican.es/tds5/dodsC/system4/System4_Seasonal_15Members.ncml",
    145 +                                      standard.vars = TRUE, dictionary = "datasets/forecasts/System4/System4_Seasonal_15Members.dic",
    146 +                                      var = "tas", members = 1,
    147 +                                      lonLim = -3.81, latLim = 43.43,
    148 +                                      season = 1, years = 1990:1999, leadMonth = 1)
    149 > madridData <- loadSeasonalForecast(dataset = "http://www.meteo.unican.es/tds5/dodsC/system4/System4_Seasonal_15Members.ncml",
    150 +                                      standard.vars = TRUE, dictionary = "datasets/forecasts/System4/System4_Seasonal_15Members.dic",
    151 +                                      var = "tas", members = 1,
    152 +                                      lonLim = -3.68, latLim = 40.40,
    153 +                                      season = 1, years = 1990:1999, leadMonth = 1)
    154 > plot(santanderData$MemberData[[1]], ty='l', col = "red")
    155 > lines(madridData$MemberData[[1]], ty='l')
    156 > title("Same data as the previous plot")
    157 
    158 }}}
    159 
    160 
    161 [[Image(timeSeriesS4_2.png)]]
    162 
    163 
    164 
    165 
    166 
    167 = loadObservations =
    168 
    169 The function `loadObservations` is intended to deal with observational datasets from weather stations stored as csv files in a standard format.
    170 In the directory ''meteoR/datasets/observations/Iberia_ECA'' there is an example dataset.
    171 
    172 {{{
    173 > list.files("./datasets/observations/Iberia_ECA")
    174 [1] "ecaIberia.nc" "Master.csv"   "pr.csv"       "tas.csv"      "tasmax.csv"   "tasmin.csv"
    175 
    176 }}}
    177 
    178 As we can see, there is a number of files with the name of the variables they store, and a ''Master.csv'' file, which contains the required metadata in order to identify each station.
    179 This is how the ''Master.csv'' file looks like:
    180 {{{     
    181 > master <- read.csv("./datasets/observations/Iberia_ECA/Master.csv")
    182 > str(master)
    183 'data.frame':      28 obs. of  5 variables:
    184 $ Id       : int  33 229 230 231 232 233 234 236 309 336 ...
    185 $ Longitude: num  1.38 -6.83 -3.68 -4.49 -4.01 ...
    186 $ Latitude : num  43.6 38.9 40.4 36.7 40.8 ...
    187 $ Altitude : int  151 185 667 7 1894 790 251 44 43 704 ...
    188 $ Metadata : Factor w/ 1 level " Data provided by the ECA&D project. Available at http://www.ecad.eu": 1 1 1 1 1 1 1 1 1 1 ...
    189 
    190 }}}
    191 
    192 The ''Master'' table contains the following fields:
    193 * `Id`: Identification code of the station. This code is used as argument by the function `loadObservations`. Note thas this field should actually be read as a character string, as internally done by the `loadObservations` function. However, in this case `read.csv` by default interpretes it as a numeric value.
    194 * `Longitude`: longitude
    195 * `Latitude`: latitude
    196 * `Altitude`: altitude
    197 * `Metadata`: other metadata associated to the dataset
    198 
    199 First of all we will plot the stations so that we can get an idea of their geographical situation and extent:
    200 
    201 {{{
    202 > library(fields)     
    203 > plot(master[ ,2:3], asp=1, pch=15, col="red")
    204 > world(add=TRUE)
    205 
    206 }}}
    207 
    208 In order to get a vector with the correct IDs as character strings instead of numeric values, we can load again the corresponding column using the argument `colClasses = "character"`
    209 
    210 {{{
    211 > stationIDs <- read.csv("./datasets/observations/Iberia_ECA/Master.csv", colClasses = "character")[ ,1]
    212 > stationIDs
    213 [1] "000033" "000229" "000230" "000231" "000232" "000233" "000234" "000236" "000309" "000336" "000414" "000416" "000420" "000421"
    214 [15] "000788" "000800" "001392" "001398" "003904" "003905" "003907" "003908" "003922" "003928" "003936" "003947" "003948" "003949"
    215 
    216 }}}
    217 
    218 We place the labels on top of each location. There are ways to avoid the overlapping of labels in order to explore the dataset, for instance in an interactive fashion by means of the `identify` function .
    219 
    220 {{{
    221 > text(master$Longitude + .2, master$Latitude + .2, stationIDs, cex=.7)
    222 
    223 }}}
    224 
    225 
    226 [[Image(stations.png)]]
    227 
    228 In this particular example we are interested in the mean surface temperature from the cities of Santander and Madrid in January. The station codes are
    229 ''"000230"'' and ''"001392"'' for Madrid and Santander respectively. We will select the period 1990-1999.
    230 
    231 {{{
    232 > stationData <- loadObservations(source.dir="./datasets/observations/Iberia_ECA/", var="tas",
    233 +           standard.vars=FALSE, stationID=c("000230","001392"), startDate="1990-01-01", endDate="1999-12-31", season = 1)
    234 > str(stationData)
    235 List of 5
    236  $ StationID   : chr [1:2] "000230" "001392"
    237  $ LatLonCoords: num [1:2, 1:2] 40.41 43.46 -3.68 -3.82
    238   ..- attr(*, "dimnames")=List of 2
    239   .. ..$ : NULL
    240   .. ..$ : chr [1:2] "Latitude" "Longitude"
    241  $ Altitude    : int [1:2] 667 64
    242  $ Dates       : POSIXlt[1:3652], format: "1990-01-01" "1990-01-02" "1990-01-03" "1990-01-04" ...
    243  $ Data        : num [1:3652, 1:2] 112 85 90 74 101 109 110 82 86 66 ...
    244   ..- attr(*, "dimnames")=List of 2
    245   .. ..$ : NULL
    246   .. ..$ : chr [1:2] "000230" "001392"
    247 > ylimits <- c(floor(min(stationData$Data)) * .1, ceiling(max(stationData$Data) * .1))
    248 > plot(stationData$Data[ ,1], ty='n', ylim = ylimits, axes=FALSE, ylab="degC", xlab="Year")
    249 > axis(1,at = seq(1,31*11,31), labels=c(1990:1999,""))
    250 > axis(2, ylim=ylimits)
    251 > abline(v=seq(1,31*10,31), lty=2, col = "grey30")
    252 > for (i in 1:ncol(stationData$Data)) {
    253 +      lines(stationData$Data[ ,i] * .1, col=i)
    254 + }
    255 > city.names = c("Madrid","Santander")
    256 > legend("topright", paste(city.names, stationData$StationID), lty=1, col=1:2, bty="n")
    257 > title(main = "Mean surface Temperature January - Observations")
    258 
    259 }}}   
    260 
    261 [[Image(timeSeriesObs2.png)]]
    262 
    263 
    264 
    265 
    266 = Creating a dataset =
    267 
    268 Climate datasets of various types (e.g. reanalysis, RCM/GCM data...) are often stored as collections of netCDF files for each particular variable and/or time slice. These files can be either locally or remotely stored. A convenient way of dealing with this kind of datasets is the use of NcML files. A NcML file is a ​XML representation of netCDF metadata. By means of NcML it is possible to create virtual datasets by modifying and aggregating other datasets, thus providing maximum flexibility and ease of access to data stored in collections of files containing data from different variables/time slices. The function `makeNcmlDataset` is intended to deal with reanalysis, forecasts and other climate data products, often consisting of collections of netCDF files of these characteristics.
    269 
    270 In this example, we have chosen several variables commonly used in statistical downscaling applications belonging to the NCEP reanalysis, and stored in different netCDF files. These variables are stored in the following directory:
    271 
    272 {{{
    273 > list.files("./datasets/reanalysis/Iberia_NCEP/")
    274 [1] "Iberia_NCEP.dic"  "NCEP_Q850.nc"     "NCEP_SLPd.nc"     "NCEP_T850.nc"     "NCEP_Z500.nc"
    275 
    276 }}}
    277 
    278 Note that the dictionary for this dataset (''Iberia_NCEP.dic'') has also been included in the directory.
    279 
    280 The function `makeNcmlDataset` is used to conveniently aggregate the required information so that the inventory/loading functions point to the NcML rather that to the netCDF files. The following call to the function wll create the NcML file in the same directory where the netCDF files are stored:
    281 
    282 {{{
    283 > makeNcmlDataset(source.dir="datasets/reanalysis/Iberia_NCEP/", ncml.file="datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml")
    284 [2013-05-20 10:00:51]
    285 NcML file "Iberia_NCEP_dataset.ncml" created from 4 files corresponding to 4 variables
    286 Use 'dataInventory(NcML file)' to obtain a description of the dataset
    287 
    288 }}}
    289 
    290 The function creates a new NcML file in the directory specified (in this case in the working directory, as no path has been specified), and gives some information about the number of files and variables conforming the dataset. In the next section is described how to find out the different variables stored in the newly created dataset and their characteristics.
    291 
    292 = dataInventory =
    293 
    294 With the aid of the `dataInventory`function we can easily retrieve all the necessary information to access and manipulate the variables stored in a dataset.  In the following example, we get a description of the NcML dataset created in the previous section, containing several variables of the NCEP reanalysis in the Iberian Peninsula.
    295 
    296 {{{
    297 > inv.iberiaNCEP <- dataInventory("datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml")
    298 # Structure of the inventory
    299 > str(inv.iberiaNCEP)
    300 List of 4
    301  $ Q   :List of 5
    302   ..$ Description: chr "Specific humidity"
    303   ..$ DataType   : chr "float"
    304   ..$ Units      : chr "kg kg**-1"
    305   ..$ TimeStep   :Class 'difftime'  atomic [1:1] 24
    306   .. .. ..- attr(*, "tzone")= chr ""
    307   .. .. ..- attr(*, "units")= chr "hours"
    308   ..$ Dimensions :List of 4
    309   .. ..$ level:List of 3
    310   .. .. ..$ Type  : chr "Pressure"
    311   .. .. ..$ Units : chr "millibar"
    312   .. .. ..$ Values: num 850
    313   .. ..$ time :List of 3
    314   .. .. ..$ Type  : chr "Time"
    315   .. .. ..$ Units : chr "days since 1950-01-01 00:00:00"
    316   .. .. ..$ Values: POSIXlt[1:16071], format: "1958-01-01" "1958-01-02" "1958-01-03" "1958-01-04" ...
    317   .. ..$ lat  :List of 3
    318   .. .. ..$ Type  : chr "Lat"
    319   .. .. ..$ Units : chr "degrees north"
    320   .. .. ..$ Values: num [1:6] 35 37.5 40 42.5 45 47.5
    321   .. ..$ lon  :List of 3
    322   .. .. ..$ Type  : chr "Lon"
    323   .. .. ..$ Units : chr "degrees east"
    324   .. .. ..$ Values: num [1:9] -15 -12.5 -10 -7.5 -5 -2.5 0 2.5 5
    325  $ SLPd:List of 5
    326   ..$ Description: chr "Mean Sea Level Pressure; Mean daily value"
    327   ..$ DataType   : chr "float"
    328   ..$ Units      : chr "Pa"
    329   ..$ TimeStep   :Class 'difftime'  atomic [1:1] 24
    330   .. .. ..- attr(*, "tzone")= chr ""
    331   .. .. ..- attr(*, "units")= chr "hours"
    332   ..$ Dimensions :List of 3
    333   .. ..$ time:List of 3
    334   .. .. ..$ Type  : chr "Time"
    335   .. .. ..$ Units : chr "days since 1950-01-01 00:00:00"
    336   .. .. ..$ Values: POSIXlt[1:16071], format: "1958-01-01" "1958-01-02" "1958-01-03" "1958-01-04" ...
    337   .. ..$ lat :List of 3
    338   .. .. ..$ Type  : chr "Lat"
    339   .. .. ..$ Units : chr "degrees north"
    340   .. .. ..$ Values: num [1:6] 35 37.5 40 42.5 45 47.5
    341   .. ..$ lon :List of 3
    342   .. .. ..$ Type  : chr "Lon"
    343   .. .. ..$ Units : chr "degrees east"
    344   .. .. ..$ Values: num [1:9] -15 -12.5 -10 -7.5 -5 -2.5 0 2.5 5
    345  $ T   :List of 5
    346   ..$ Description: chr "Temperature"
    347   ..$ DataType   : chr "float"
    348   ..$ Units      : chr "K"
    349   ..$ TimeStep   :Class 'difftime'  atomic [1:1] 24
    350   .. .. ..- attr(*, "tzone")= chr ""
    351   .. .. ..- attr(*, "units")= chr "hours"
    352   ..$ Dimensions :List of 4
    353   .. ..$ level:List of 3
    354   .. .. ..$ Type  : chr "Pressure"
    355   .. .. ..$ Units : chr "millibar"
    356   .. .. ..$ Values: num 850
    357   .. ..$ time :List of 3
    358   .. .. ..$ Type  : chr "Time"
    359   .. .. ..$ Units : chr "days since 1950-01-01 00:00:00"
    360   .. .. ..$ Values: POSIXlt[1:16071], format: "1958-01-01" "1958-01-02" "1958-01-03" "1958-01-04" ...
    361   .. ..$ lat  :List of 3
    362   .. .. ..$ Type  : chr "Lat"
    363   .. .. ..$ Units : chr "degrees north"
    364   .. .. ..$ Values: num [1:6] 35 37.5 40 42.5 45 47.5
    365   .. ..$ lon  :List of 3
    366   .. .. ..$ Type  : chr "Lon"
    367   .. .. ..$ Units : chr "degrees east"
    368   .. .. ..$ Values: num [1:9] -15 -12.5 -10 -7.5 -5 -2.5 0 2.5 5
    369  $ Z   :List of 5
    370   ..$ Description: chr "Geopotential"
    371   ..$ DataType   : chr "float"
    372   ..$ Units      : chr "m**2 s**-2"
    373   ..$ TimeStep   :Class 'difftime'  atomic [1:1] 24
    374   .. .. ..- attr(*, "tzone")= chr ""
    375   .. .. ..- attr(*, "units")= chr "hours"
    376   ..$ Dimensions :List of 4
    377   .. ..$ level:List of 3
    378   .. .. ..$ Type  : chr "Pressure"
    379   .. .. ..$ Units : chr "millibar"
    380   .. .. ..$ Values: num 850
    381   .. ..$ time :List of 3
    382   .. .. ..$ Type  : chr "Time"
    383   .. .. ..$ Units : chr "days since 1950-01-01 00:00:00"
    384   .. .. ..$ Values: POSIXlt[1:16071], format: "1958-01-01" "1958-01-02" "1958-01-03" "1958-01-04" ...
    385   .. ..$ lat  :List of 3
    386   .. .. ..$ Type  : chr "Lat"
    387   .. .. ..$ Units : chr "degrees north"
    388   .. .. ..$ Values: num [1:6] 35 37.5 40 42.5 45 47.5
    389   .. ..$ lon  :List of 3
    390   .. .. ..$ Type  : chr "Lon"
    391   .. .. ..$ Units : chr "degrees east"
    392   .. .. ..$ Values: num [1:9] -15 -12.5 -10 -7.5 -5 -2.5 0 2.5 5
    393 
    394 }}}
    395 
    396 
    397 As we can see, the inventory consists of a list of four elements, which are the four variables stored in the dataset:
    398 
    399 {{{
    400 > names(inv.iberiaNCEP)
    401 [1] "Q"    "SLPd" "T"    "Z" 
    402 
    403 }}}
    404 
    405 
    406 This is how we can check the spatial domain of the dataset:
    407 
    408 {{{
    409 > lats = inv.iberiaNCEP$SLPd$Dimensions$lat$Values
    410 > lons = inv.iberiaNCEP$SLPd$Dimensions$lon$Values
    411 > plot(expand.grid(lons,lats), asp=1)
    412 # 'world' from library 'fields'
    413 > world(add=TRUE)
    414 
    415 }}}
    416 
    417 [[Image(iberiaNCEPextent.png)]]
    418 
    419 = loadGCM
    420 
    421 Once the NcML dataset is created and we get an idea of the nature of the variables stored, the `loadGCM` function is used to retrieve the variables desired at selected dimensional slices. Although the name of the function may result somewhat misleading, the function is intended for loading many kinds of gridded datasets, and not only GCM data, including reanalysis, RCM data and observational gridded datasets, for instance. In this particular example, we will load the temperature data from the NCEP reanalysis in the Iberian Peninsula, provided in the example datasets of the meteoR package.
    422 
    423 We have a look again to the description of the variable temperature, as provided by the `dataInventory`:
    424 
    425 {{{
    426 > str(inv.iberiaNCEP$T)
    427 List of 5
    428  $ Description: chr "Temperature"
    429  $ DataType   : chr "float"
    430  $ Units      : chr "K"
    431  $ TimeStep   :Class 'difftime'  atomic [1:1] 24
    432   .. ..- attr(*, "tzone")= chr ""
    433   .. ..- attr(*, "units")= chr "hours"
    434  $ Dimensions :List of 4
    435   ..$ level:List of 3
    436   .. ..$ Type  : chr "Pressure"
    437   .. ..$ Units : chr "millibar"
    438   .. ..$ Values: num 850
    439   ..$ time :List of 3
    440   .. ..$ Type  : chr "Time"
    441   .. ..$ Units : chr "days since 1950-01-01 00:00:00"
    442   .. ..$ Values: POSIXlt[1:16071], format: "1958-01-01" "1958-01-02" "1958-01-03" "1958-01-04" ...
    443   ..$ lat  :List of 3
    444   .. ..$ Type  : chr "Lat"
    445   .. ..$ Units : chr "degrees north"
    446   .. ..$ Values: num [1:6] 35 37.5 40 42.5 45 47.5
    447   ..$ lon  :List of 3
    448   .. ..$ Type  : chr "Lon"
    449   .. ..$ Units : chr "degrees east"
    450   .. ..$ Values: num [1:9] -15 -12.5 -10 -7.5 -5 -2.5 0 2.5 5
    451 
    452 }}}
    453 
    454 As we can see, the variable T has vertical levels. In this case, the only level available is at 850 mb. The variable is daily, as we can see in the `TimeStep` element of the list, and the original units are Kelvin.
    455 
    456 There are several options for spatial selection using the `loadGCM` function, as in the case of `loadSeasonalForecast`. For instance, if we want the whole domain of the dataset, there is no need for specifying the `lonLim` and `latLim` arguments. Alternatively, it is possible to select smaller rectangular domains or single points. In the next example, we will load T850 for January in a similar domain than previously with the System4 dataset, centered on the Iberian Peninsula, encompassing the period 1990-1999.
    457 
    458 
    459 {{{
    460 > t850.ncep.iberia <- loadGCM(dataset = "./datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml", standard.vars=TRUE,
    461 +        var="ta", lonLim = c(-10,5), latLim = c(35,45), level=850, season=1, years=1990:1999)
    462 > str(t850.ncep.iberia)
    463 List of 5
    464  $ VarName     : chr "ta"
    465  $ Level       : num 850
    466  $ Dates       :List of 2
    467   ..$ Start: POSIXlt[1:310], format: "1990-01-01" "1990-01-02" "1990-01-03" "1990-01-04" ...
    468   ..$ End  : POSIXlt[1:310], format: "1990-01-01" "1990-01-02" "1990-01-03" "1990-01-04" ...
    469  $ LatLonCoords: num [1:35, 1:2] 35 37.5 40 42.5 45 35 37.5 40 42.5 45 ...
    470   ..- attr(*, "dimnames")=List of 2
    471   .. ..$ : NULL
    472   .. ..$ : chr [1:2] "lat" "lon"
    473  $ Data        : num [1:310, 1:35] 5.95 5.55 0.35 2.05 5.35 ...
    474 
    475 }}}
    476 
    477 Note that in this particular case, we are loading standard variables, as defined in the vocabulary (by setting the argument `standard.vars = TRUE`), but we did not specify a path to the dictionary (default to `NULL`). By default, the function searches the dictionary in the same directory where the ''ncml'' file has been created, assuming that this is a file with extension ''.dic'' and the same name as the ''ncml''.
    478 
    479 
    480 The function `loadGCM` allows the retrieval of rectangular windows as well as single point locations. In the last case, the function finds the closest grid point to the given coordinates, as previously described in the example of `loadSeasonalForecast`. In the following lines we extract the time series of the Spanish cities used in the previous examples:
    481 
    482 {{{
    483 > t850.ncep.santander <- loadGCM(dataset = "./datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml", standard.vars=TRUE,
    484 +                            var="ta", lonLim = -3.81, latLim = 43.43, level=850, season=1, years=1990:1999)
    485 > t850.ncep.madrid <- loadGCM(dataset = "./datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml", standard.vars=TRUE,
    486 +                               var="ta", lonLim = -3.68, latLim = 40.40, level=850, season=1, years=1990:1999)
    487 >
    488 > ylimits <- c(floor(min(t850.ncep.madrid$Data)), ceiling(max(t850.ncep.santander$Data)))
    489 > plot(t850.ncep.santander$Data, ty='n', axes=FALSE, ylab="degC", xlab="Year", ylim = ylimits)
    490 > axis(1,at = seq(1,31*11,31), labels=c(1990:1999,""))
    491 > axis(2, ylim=ylimits)
    492 > abline(v=seq(1,31*10,31), lty=2)
    493 > lines(t850.ncep.santander$Data, ty='l')
    494 > lines(t850.ncep.madrid$Data, col = "red")
    495 > legend("bottomleft", city.names, lty=1, col=1:4)
    496 > title(main = "Mean Temperature January 850mb")
    497 > mtext("NCEP reanalysis")
    498 
    499 }}}
    500 
    501 
    502 [[Image(timeSeriesNCEP.png)]]
     26spplot(sgdf, zcol = 1:4, scales=list(draw = TRUE), col.regions = rev(terrain.colors(50)), at = seq(0, ceiling(max(sgdf@data)),10), sp.layout = list(l1))