spnc: Spatial and NetCDF
spnc provides a uniform interface to NetCDF files for accessing data in the
sp framework of classes. Using Unidata's netcdf_4 library permits access to both locally stored data (as .nc files) as well as network served data (as OpeNDAP/DODS/Thredds/etc/).
Under the hood we use David Pierce's ncdf4 R package to interface with the Unidata netcdf_4 library. We try to keep the underlying resources opaque to the user so that we may switch to another interface to Unidata's netcdf_4 library as the need arises.
R packages... may require substantial effort to install on some platforms
The repository is set up for installation using the devtools package.
library(devtools) # for your own library install_github("btupper/spnc") # or for system wide installation (requires admin-level permissions) with_lib(.Library, install_github("btupper/spnc"))
If you are developing new data source readers and have the repo on your local platform then use the following.
library(devtools) pkg_path <- '/path/to/spnc' install(pkg_path)
If you edit the code (with documentation as needed!) then you can easily re-install - sort of on the fly.
library(spnc) x <- SPNC(some_resource) # edit, document and then redo the document/install process library(devtools) pkg_path <- '/path/to/spnc' document(pkg_path) install(pkg_path) # you may need to reinstantiate your reference *if* the definition or methods have changed in your edits x <- SPNC(some_resource)
We try to follow the KISS principle by minimizing the exposure of details. See the wiki for examples.
SPNC(ncdf_resource_or_filename, bb)create an instance of SPNCRefClass or subclass with the provided bounding box.
bb is a 4 element vector that follows the pattern of the
Extent class in the
sp package: [xmin, xmax, ymin, ymax]. If you set the bounding box when you instantiate the SPNC object, then it is used as the default for subsequent calls to the methods shown below, unless you specifically override the value.
open()open the NetCDF connection
close()close the NetCDF connection
flavor()retrieve the flavor of the NetCDF data (e.g. [source='OISST',type='raster', local = TRUE/FALSE/NA])
get_raster(bb, t)retrieve a n-d array as SpatialRaster within the bounding box for the given times
get_path(bb, t)retrieve a path SpatialLinesDataFrame within the bounding box for the given times (not sure what this is just yet!)
get_points(bb, t)retrieve set of points as SpatialPointsDataFrame within the bounding box for the given times
Known sources of data
The Ocean Color group provides Level3 Standard Mapped Images (L3SMI) for a variety of products. See more in the wiki on how to programmatically search their database. Tested on MODISA CHL and SST. Note, a companion package, obpgcrawler, is available to programmatically find Ocean Color data.
JPL's PO.DAAC provides this data as an NCML (filename.ncml) - an XML file that the NetCDF-4 library happily reads as if it were a NetCDF file. THis provides easy access to data at different times, so not only does one collect data from geographic subregions but for different times. See the wiki for details.
Currently this will access (a) locally downloaded files and (b) OpeNDAP files for a single time (dailies). Access to the aggregate NCML OpeNDAP resources is still under development. Note that longitude is stored in the 0-360 format which requires some juggling internally.
Note that you must specify bounding boxes in [0,360] longitude format. OISST is stored with longitudes referenced an the range [0,360]. Use
spnc::to360BB() to convert longitudes from [-180,180] to [0,360].
Bwlow we access the same data two ways...
- (a) download and unpack this file
- (b) connect via OpeNDAP to this resource
... and visually compare the two (which are the same!)
library(raster) library(spnc) # note the transform from [-180,180] to [0,360] longitude BB <- spnc::to360BB(c(-72,-63,39,46)) localFile <- '/Users/ben/Downloads/amsr-avhrr-v2.20040101.nc' X <- SPNC(localFile) x <- X$get_raster(what = 'sst', bb = BB) opendapFile <- 'http://www.ncdc.noaa.gov/thredds/dodsC/oisst/NetCDF/AVHRR-AMSR/2004/AVHRR-AMSR/amsr-avhrr-v2.20040101.nc' Y <- SPNC(opendapFile) y <- Y$get_raster(what = 'sst', bb = BB) spplot(stack(x,y)) # for NCML mytimes <- seq(from = as.POSIXct("1985-01-10", tz = 'UTC'), to = as.POSIXct("1997-01-01", tz = 'UTC'), by = 'year') OISST_file = 'http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR_agg' OI <- SPNC(OISST_file, bb = BB) oi <- OI$get_raster(what = 'sst', time = mytimes) spplot(oi)
See more at the wiki
library(raster) library(spnc) # note the transform from [-180,180] to [0,360] longitude BB <- spnc::to360BB(c(-72,-63,39,46)) ff <- 'http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/cpc_us_precip/RT/precip.V1.0.2015.nc' X <- SPNC(ff) r <- X$get_raster(bb = BB, time = 1:5) spplot(r)
NARR is provided for the continental US plus a generous extension. Data are stored with project coordinates so there are some signifiant difference between this and other classes that inherit from
SPNCRefClass - see this for more info. Perhap the most significant difference is that the data are stored in a
raster::RasterLayer form rather than as ncdf. You might consider not using SPNC to access NARR datasets. The 'best' utility of
NARRRefClass is that requests for points or subset rasters can be done using long/lat values rather than the projected coordinates.
library(raster) library(spnc) ff <- '/Users/Shared/data/ecocast/wx/air.sfc.2015.nc' bb <- c(-72,-63,39,46) X <- SPNC(ff) R <- X$get_raster() spplot(R) r <- X$get_raster(bb = bb) spplot(r)
Blended Sea WindsNOAA/NCEI
Blended Sea Winds provides global daily coverage of
w winds on a 0.25 degree going back the 1980s.
library(raster) library(rasterVis) library(spnc) uri <-'https://www.ncdc.noaa.gov/thredds/dodsC/oceanwindsdly' BB <- c(-72,-63,39,46) BB360 <- spnc::to360BB(BB) X <- SPNC(uri, bb = BB360) layer = as.POSIXct("2015-08-15", tz = 'UTC') u <- X$get_raster(what = 'u', layer = layer) v <- X$get_raster(what = 'v', layer = layer) R <- raster::stack(u,v) R <- raster::rotate(R) names(R) <- c("u", "v") P <- rasterVis::vectorplot(R, isField = 'dXY') png(); print(P); dev.off()
snow <- SPNC('http://www.ncdc.noaa.gov/thredds/dodsC/cdr/snowcover/nhsce_v01r01_19661004_latest.nc')
Under development. Provides mask and areal extent using a modified grid.