Recipe: R program as OPeNDAP client: Difference between revisions
Jump to navigation
Jump to search
Initial add |
mNo edit summary |
||
| Line 9: | Line 9: | ||
== Solution == | == Solution == | ||
The following code has been tested with R version 3.1.2. | |||
The metadata for the collection <code>inst01hr_3d_T_Cv</code> is available at | The metadata for the collection <code>inst01hr_3d_T_Cv</code> is available at | ||
| Line 26: | Line 26: | ||
library(ncdf4) | library(ncdf4) | ||
library(rworldmap) | library(rworldmap) | ||
# Url | # Url | ||
url <- "http://opendap.nccs.nasa.gov:80/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv" | url <- "http://opendap.nccs.nasa.gov:80/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv" | ||
| Line 45: | Line 42: | ||
t <- ncvar_get(nc,"t",start=offset,count=count) | t <- ncvar_get(nc,"t",start=offset,count=count) | ||
# Statistics | # Statistics | ||
summary(t) | print(summary(t)) | ||
# Plot Data | # Plot Data | ||
mapGriddedData(t[1:im,1:jm,71],xlim=c(-130.,-65.),ylim=c(25.,50.)) | mapGriddedData(t[1:im,1:jm,71],xlim=c(-130.,-65.),ylim=c(25.,50.)) | ||
| Line 54: | Line 49: | ||
==== Output ==== | ==== Output ==== | ||
Running this | Running this R script | ||
> source("g5nr.R") | |||
we get the text output | we get the text output | ||
Min. 1st Qu. Median Mean 3rd Qu. Max. | |||
190.7 221.2 248.1 248.7 275.3 305.8 | |||
and the plot | and the plot | ||
[[Image:SurfaceAirTempSubsetIDL.png|500px]] | <!-- [[Image:SurfaceAirTempSubsetIDL.png|500px]] --> | ||
== Discussion == | == Discussion == | ||
Latest revision as of 10:14, 9 March 2015
Back to G5NR Data Access Guide.
Problem
By accessing the collection inst01hr_3d_T_Cv via the OPeNDAP server
http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km
we want to read the surface temperature data inside the box bound by latitudes 25oN, 50oN and longitudes -130oW, -65oW for 2006/Sep/18, 9z, compute its min/max and plot the data at the surface.
Solution
The following code has been tested with R version 3.1.2.
The metadata for the collection inst01hr_3d_T_Cv is available at
http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv.info
Code
This code accesses the collection inst01hr_3d_T_Cv from the OPeNDAP server and reads a subset of the temperature data (all levels inside the bounding box specified above) and computes its max/min. It then plots the data at the surface (level=72).
NOTE:
- Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.
- While in the downloaded file, the temperature variable appears in the uppercase (T), on the OPeNDAP server, this variable is in lowercase.
- Via the OPeNDAP URL, we now have access to all times for which data exists. The hourly
instfiles are available starting at 2005/15/15, 2200z. Our desired time, 2006/09/18, 0900z is then the 11772th step.
# Load Libraries
library(ncdf4)
library(rworldmap)
# Url
url <- "http://opendap.nccs.nasa.gov:80/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv"
nc <- nc_open(url)
# Bounding Box
imin <- round((-130.+180.)/0.5)-1
imax <- round(( -64.+180.)/0.5)-1
jmin <- round(( 25 + 90.)/0.5) - 1
jmax <- round(( 50 + 90.)/0.5) - 1
im <- imax-imin+1
jm <- jmax-jmin+1
lm <- 72
# Start and Count vectors
offset <- c(imin, jmin, 1, 11771)
count <- c(im, jm, lm, 1)
t <- ncvar_get(nc,"t",start=offset,count=count)
# Statistics
print(summary(t))
# Plot Data
mapGriddedData(t[1:im,1:jm,71],xlim=c(-130.,-65.),ylim=c(25.,50.))
Output
Running this R script > source("g5nr.R")
we get the text output
Min. 1st Qu. Median Mean 3rd Qu. Max. 190.7 221.2 248.1 248.7 275.3 305.8
and the plot
Discussion
See Also
- File Spec: File:G5NR-Ganymed-7km FileSpec-ON6-V1.0.pdf
- Recipe: Fortran program to read data from downloaded file