wiki:udg/ecoms/RPackage/examples

Version 13 (modified by juaco, 9 years ago) (diff)

--

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.

> openDAP.query <- loadSystem4(dataset = "http://www.meteo.unican.es/tds5/dodsC/system4/System4_Seasonal_15Members.ncml", 
+                             var = "tas", members = 1,
+                             lonLim = c(-10,5), latLim = c(35,45),
+                             season = 1, years = 1990:1999, leadMonth = 1)

Data are now loaded into the R session:

> str(openDAP.query)
List of 7
 $ VarName      : chr "Mean_temperature_at_2_metres"
 $ VarUnits     : chr "degC"
 $ TimeStep     :Class 'difftime'  atomic [1:1] 1
  .. ..- attr(*, "tzone")= chr ""
  .. ..- attr(*, "units")= chr "days"
 $ MemberData   :List of 1
  ..$ : num [1:310, 1:280] 13.3 13.9 12.5 13 13 ...
 $ LatLonCoords : num [1:280, 1:2] 45 44.2 43.5 42.7 42 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "lat" "lon"
 $ RunDates     : POSIXlt[1:310], format: "1989-12-01" "1989-12-01" "1989-12-01" ...
 $ ForecastDates:List of 2
  ..$ Start: POSIXlt[1:310], format: "1990-01-01" "1990-01-02" "1990-01-03" ...
  ..$ End  : POSIXct[1:310], format: "1990-01-02" "1990-01-03" "1990-01-04" ...

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:

> city.names <- c("Sevilla", "Madrid", "Santander", "Zaragoza")
> locations <- matrix(c(-5.9, 37.4167, -3.68, 40.4, -3.817, 43.43, -0.8167, 41.667), ncol=2, byrow = TRUE)
> dimnames(locations) <- list(city.names, c("lon","lat"))
> print(locations)
             lon     lat
Sevilla   -5.9000 37.4167
Madrid    -3.6800 40.4000
Santander -3.8170 43.4300
Zaragoza  -0.8167 41.6670

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:

> # Representation of the mean temperature of the period
> mean.field <- colMeans(openDAP.query$MemberData[[1]])
> lat <- openDAP.query$LatLonCoords[ ,1]
> lon <- openDAP.query$LatLonCoords[ ,2]
> # Requires "akima::interp" for regular grid interpolation 
> library(akima)
> filled.contour(interp(lon, lat, mean.field), asp=1,
+               plot.title = title (main = "Mean surface T January 1990-99",ylab = "latitude", xlab = "longitude"),
+               key.title = title(main = "degC"),
+               key.axes = axis(4, seq(-1,16,1)),
+               plot.axes = {points(locations, pch=15); axis(1); axis(2); text(locations, cex=.8, pos=3, city.names)},
+               color.palette = topo.colors)

No image "TmeanJan.png" attached to udg/ecoms/RPackage/examples

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.

> # Selection of point locations. Requires "sp::spDistsN1" to compute geographic distances
> library(sp)
> index <- rep(NA, nrow(locations))
> for (i in 1:length(index)) {
+      index[i] <- which.min(spDistsN1(openDAP.query$LatLonCoords[ ,2:1], locations[i, ]))      
+ }
> locations.data <- openDAP.query$MemberData[[1]][ ,index]
> colnames(locations.data) <- city.names
> str(locations.data)
 num [1:310, 1:4] 7.71 8.78 10.77 10.88 11.55 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:4] "Sevilla" "Madrid" "Santander" "Zaragoza"

The object locations.data is a matrix in which time series are arranged in columns for each of the four locations selected.

> ylimits <- c(floor(min(locations.data)), ceiling(max(locations.data)))
> plot(locations.data[ ,1], ty='n', ylim = ylimits, axes=FALSE, ylab="degC", xlab="Year")
> axis(1,at = seq(1,31*11,31), labels=c(1990:1999,""))
> axis(2, ylim=ylimits)
> abline(v=seq(1,31*10,31), lty=2)
> for (i in 1:ncol(locations.data)) {
+      lines(locations.data[ ,i], col=i)
+ }
> legend("bottomleft", city.names, lty=1, col=1:4)
> title(main = "Mean surface Temperature January")

No image "timeSeries.png" attached to udg/ecoms/RPackage/examples