In bgctw/ncdfTools: Handling NetCDF Files

# do not execute on CRAN: 
# https://stackoverflow.com/questions/28961431/computationally-heavy-r-vignettes
is_check <- ("CheckExEnv" %in% search()) || any(c("_R_CHECK_TIMINGS_",
             "_R_CHECK_LICENSE_") %in% names(Sys.getenv()))
# knitr::opts_chunk$set(eval = !is_check)

#rmarkdown::render("vignettes/useCase.Rmd")
knitr::opts_knit$set(root.dir = '..')
knitr::opts_chunk$set(
    #, fig.align = "center"
    #, fig.width = 3.27, fig.height = 2.5, dev.args = list(pointsize = 10)
    #,cache = TRUE
    #, fig.width = 4.3, fig.height = 3.2, dev.args = list(pointsize = 10)
    #, fig.width = 6.3, fig.height = 6.2, dev.args = list(pointsize = 10)
    # works with html but causes problems with latex
    #,out.extra = 'style = "display:block; margin: auto"' 
    )
knitr::knit_hooks$set(spar = function(before, options, envir) {
    if (before) {
        par(las = 1 )                   #also y axis labels horizontal
        par(mar = c(2.0,3.3,0,0) + 0.3 )  #margins
        par(tck = 0.02 )                          #axe-tick length inside plots             
        par(mgp = c(1.1,0.2,0) )  #positioning of axis title, axis labels, axis
     }
})

#themeTw <- theme_bw(base_size = 10) + theme(axis.title = element_text(size = 9))
#bgiDir <- "~/bgi"

Typical workflow

Creating a new point-nc-File

require(ncdfTools)
nRec <- 30L
times <- seq(ISOdatetime(2010,1,1,0,0,0, tz = "UTC")
             , by = "30 min", length.out = nRec)
dsTest <- data.frame(
  time = times
  , soilResp = rnorm(nRec)
  , sdSoilResp = 0.1
)
dsTest$sdSoilResp[2:4] <- NA
tmpDir <- tempdir()
#tmpDir <- "tmp"
ncFileDs <- createStdNcdfFile(
  data = dsTest
  , file.name = file.path(tmpDir,"dsSoilResp.nc")
  , lat.values = 53.0
  , long.values = 13
  , time.values = times
  , units = rep('g/m^2',2)
  )
# single variable: can supply data as vector
ncFile <- createStdNcdfFile(
  "soilResp"
  , file.name = file.path(tmpDir,"soilResp.1.1.30.nc")
  , lat.values = 53.0
  , long.values = 13
  , time.values = dsTest$time
  , units = 'g/m2'
  , data = dsTest$soilResp
  )
fileName <- ncFile

Getting information on file structure

Dimensions

infoNcdfDims(ncFile)

Variables

#infoNcdfVars(ncFile, order.var = 'id')
infoNcdfVars(ncFileDs, order.var = 'id')

global attributes

as.data.frame(infoNcdfAtts(ncFile))

and variable attributes

as.data.frame(infoNcdfAtts(ncFile, "soilResp"))

Coordinate values

readCoordinatesNcdf(ncFile)

Try matching the main variable (not dimensions, having most dimenstion, matching filename)

readNcdfVarName(ncFile)

typing-sparse reading of a single variable from ncdf-File

ans <- readNcdf(ncFile) 
all(as.vector(ans) - dsTest$soilResp < 1e-10 )

checking file structure compliance to CF/COARDS/MDI ncdf conventions

checkNcdfFile(ncFile)

Modifying NetCDF files in place

appending string to the history

modifyNcdfAppendHistory( ncFile, "another history entry")
cat(infoNcdfAtts(ncFile)["history","value"])

see also modifyNcdfDefAtts, modifyNcdfCopyAtts, modifyNcdfCopyMetadata, modifyNcdfCopyVar.

Conversion between R and NetCDF

Dates in netcdf are stored as "days since 1582-10-15 00:00". Hence, POSIXt values need to be converted to this number.

ncCon <- RNetCDF::open.nc(fileName)
timesNCOrig <- RNetCDF::var.get.nc(ncCon, 'time')
RNetCDF::close.nc(ncCon)
head(timesR <- readNcdfTime(fileName))
head(POSIX2ncdf(timesR)) - POSIX2ncdf(timesR[1])
all( (POSIX2ncdf(timesR) - timesNCOrig) == 0)

CDO commands

time aggregation. In the example take the mean over each 2 time steps.

tmpDir <- tempdir()
ans <- aggregateNcdf(ncFile, path.out = tmpDir, period = 2)
#infoNcdfDims(ncFile)
#infoNcdfDims(ans)

Adding and deleting variables

Here, functions from RNetCDF package are used. First the variable, needs to be created, this is done easiest by copying an existing one, so that dimensions and attributes are taken care of.

require(RNetCDF)
varName <- "soilRespSquared"
# work on tempfile to not change the original file
tmpFile <- tempfile(); tmp <- file.copy(fileName, tmpFile)
ncCon <- open.nc(tmpFile, write = TRUE)
tryCatch({
  resp <- var.get.nc(ncCon, 'soilResp')
  ans <- modifyNcdfCopyVar(
    ncCon, var.id.orig = "soilResp", var.id.copy = varName)
  #TODO change unit
  respSquared <- resp^2
  var.put.nc(ncCon, 'soilRespSquared', respSquared, count = c(1,1,length(resp)))
}, finally = {
 close.nc(ncCon)
})
infoNcdfVars(tmpFile)  # note line with variable soilRespSquared
# delete variable again (uses CDO command)
tmp <- modifyNcdfDeleteVar(tmpFile, varNames = varName)
infoNcdfVars(tmpFile)  # note that soilRespSquared disappeared

Reading and writing data.frames

First write a data.frame to file. Note that colum 'time' by default is written to the time dimension.

tmpFile <- tempfile()
str(dsTest)
fName <- createStdNcdfFile(
  data = dsTest
  , file.name = tmpFile
  , lat.values = 53.0
  , long.values = 13
  , units = rep('g/m2',2)
  )
infoNcdfVars(fName, dimvars = TRUE)

Use RNetCDF to read data into a general list.

require(RNetCDF)
ncCon <- open.nc(fName)
tryCatch({
  ans <- read.nc(ncCon)
}, finally = {
 close.nc(ncCon)
})
str(ans)

Use ncdfTools to read single dimensional data into a data.frame. The default single dimension is time and is by default converted to POSIXct.

ans <- ds1 <- readNcdfDataframe(fName)
head(ans)

Use ncdfTools to update a single dimensional data by a data.frame. The name of the dimension and the corresponding column in the data.frame must match.

dsMod <- ds1[5:10, c("time","sdSoilResp")]
dsMod$sdSoilResp <- 0.2
updateNcdfDataframe(dsMod, fName)
head(readNcdfDataframe(fName))

bgctw/ncdfTools documentation built on Jan. 29, 2020, 1:16 p.m.