Changes between Version 13 and Version 14 of udg/ecoms/RPackage/examples


Ignore:
Timestamp:
Apr 19, 2013 1:32:55 PM (9 years ago)
Author:
juaco
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • udg/ecoms/RPackage/examples

    v13 v14  
    1 In the next lines we describe an illustrative example of the `loadSystem4` function. 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, running this example took approximately 21 seconds.
     1In the next lines we describe an illustrative example of the `loadSystem4` function. 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:
    22
    33{{{
    … …  
    88}}}
    99
    10 Data are now loaded into the R session:
     10Data are now ready for analysis into our R session:
    1111
    1212{{{
    … …  
    3030}}}
    3131
    32 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. As an example, we will display time series of the requested dataset at four grid points coincident with the locations of four Spanish cities:
     32A 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 of the requested dataset at four grid points coincident with the locations of four 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.
    3333
    3434{{{
     35> # Calculation of the mean values of the period for each grid cell:
     36> mean.field <- colMeans(openDAP.query$MemberData[[1]])
     37> # Creation of a matrix with selected point locations:
    3538> city.names <- c("Sevilla", "Madrid", "Santander", "Zaragoza")
    3639> locations <- matrix(c(-5.9, 37.4167, -3.68, 40.4, -3.817, 43.43, -0.8167, 41.667), ncol=2, byrow = TRUE)
    … …  
    4447}}}
    4548
    46 In the following lines of code we calculate the mean temperature field and plot it, In addition, we also add to the map the point locations of the selected cities for which the time series will be represented:
     49Note 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 `akima` provides a extremely fast interpolation algorithm by means of the `interp` function.
     50
     51{{{
     52> lat <- openDAP.query$LatLonCoords[ ,1]
     53> lon <- openDAP.query$LatLonCoords[ ,2]
     54# Requires "akima::interp" for regular grid (bilinear) interpolation
     55> library(akima)
     56> grid.075 <- interp(lon, lat, mean.field, xo = seq(min(lon), max(lon), .75), yo = seq(min(lat), max(lat), .75))
     57}}}
     58
     59After 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`...)
     60
     61{{{
     62> str(grid.075)
     63List of 3
     64 $ x: num [1:19] -9.75 -9 -8.25 -7.5 -6.75 ...
     65 $ y: num [1:14] 35.2 36 36.7 37.5 38.2 ...
     66 $ z: num [1:19, 1:14] 15.1 15.1 15 15 14.3 ...
     67}}}
     68
     69In 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 `fields` provides many useful tools for spatial data handling and representation, including a world map that can be easily incorporated in our plots.
    4770
    4871{{{
    4972> # Representation of the mean temperature of the period
    50 > mean.field <- colMeans(openDAP.query$MemberData[[1]])
    51 > lat <- openDAP.query$LatLonCoords[ ,1]
    52 > lon <- openDAP.query$LatLonCoords[ ,2]
    53 > # Requires "akima::interp" for regular grid interpolation
    54 > library(akima)
    55 > filled.contour(interp(lon, lat, mean.field), asp=1,
    56 +               plot.title = title (main = "Mean surface T January 1990-99",ylab = "latitude", xlab = "longitude"),
    57 +               key.title = title(main = "degC"),
    58 +               key.axes = axis(4, seq(-1,16,1)),
    59 +               plot.axes = {points(locations, pch=15); axis(1); axis(2); text(locations, cex=.8, pos=3, city.names)},
    60 +               color.palette = topo.colors)
     73> library(fields)
     74> image.plot(grid.075, asp=1,ylab = "latitude", xlab = "longitude", col=topo.colors(36),
     75+           main = "Mean surface T January 1990-99", legend.args=list(text="degC", side=3, line=1))
     76# Adds selected locations to the plot and puts labels:
     77> points(locations, pch=15)
     78> text(locations, pos=3, city.names)
     79# Adds the world map to the current plot:
     80> world(add=TRUE)
    6181}}}
    6282
    6383[[Image(TmeanJan.png)]]
    6484
    65 Next, we plot the time series for the selected locations. To this aim, we calculate the nearest grid points to the specified locations. Note that for this task we have chosen the function `spDistsN1` from R package `sp`.
     85Next, 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 `loadSystem4` 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`):
    6686
    6787{{{
    68 > # Selection of point locations. Requires "sp::spDistsN1" to compute geographic distances
    69 > library(sp)
    70 > index <- rep(NA, nrow(locations))
    71 > for (i in 1:length(index)) {
    72 +      index[i] <- which.min(spDistsN1(openDAP.query$LatLonCoords[ ,2:1], locations[i, ]))     
     88> # Creation of a Euclidean distance matrix among all pairs of selected locations and grid points
     89> dist.matrix <- rdist(openDAP.query$LatLonCoords[ ,2:1], locations)
     90> # Positions of the nearest grid points
     91> index <- rep(NA,ncol(dist.matrix))
     92> for (i in 1:ncol(dist.matrix)) {
     93+      index[i] <- which.min(dist.matrix[ ,i])
    7394+ }
     95> # index contains the column positions in the `MemberData` matrices
    7496> locations.data <- openDAP.query$MemberData[[1]][ ,index]
    7597> colnames(locations.data) <- city.names