Cleaned up file handling (file.access for windows, ncdf4 units, nc fi…

…le handling, example code, ?)

DESCRIPTION

@@ -1,11 +1,11 @@
 Package: REddyProc
 Title: Data processing and plotting utilities of (half-)hourly eddy-covariance measurements
-Version: 0.40
+Version: 0.41
 Date: 2013-06-24
 Author: (Department for Biogeochemical Integration at MPI-BGC, Jena, Germany)
 Maintainer: Antje M. Moffat <[email protected]>
 Description: R eddy processing package with MDS gap filling algorithm, adopted after PV-Wave source code from Markus Reichstein.
 License: GPL-2
 LazyLoad: yes
 Depends: R (>= 2.13), methods
-Suggests: RNetCDF, testthat, inlinedocs
+Suggests: ncdf4, RNetCDF, testthat, inlinedocs

R/DataFunctions.R

@@ -241,7 +241,7 @@ fExpandToFullYear <- function(
 
   if( sum(!is.na(DataYear.V.n)) == 0 )
   {
-    ExpData.F.n <- data.frame(cbind(DateTime=FullYear.V.p, Data=rep(NA, length(FullYear.V.p))))
+    ExpData.F.n <- data.frame(cbind(DateTime=FullYear.V.p, Data=rep(NA_real_, length(FullYear.V.p))))
     warning(CallFunction.s, ':::fExpandToFullYear::: Variable \'', attr(Data.V.n,'varnames'), '\' contains no data for year ', Year.i, '!')
   } else if (length(TimeYear.V.p != length(FullYear.V.p))) {
     ExpData.F.n <- merge(cbind(DateTime=FullYear.V.p), cbind(DateTime=TimeYear.V.p, Data=DataYear.V.n), by='DateTime', all=T, sort=T)

R/Eddy.R

@@ -222,8 +222,10 @@ attr(sEddyProc.example,'ex') <- function(){
     EddyProc.C$sPlotDailySums('NEE_f','NEE_fsd')
     EddyProc.C$sPlotDiurnalCycle('NEE_f')
 
-    #+++ Plot other individual years/months to screen (of current R graphics device)
-    EddyProc.C$sPlotDailySumsY('NEE_f','NEE_fsd', Year.s=1998)
+    #+++ Plot individual years/months to screen (of current R graphics device)
+    EddyProc.C$sPlotHHFluxesY('NEE_f', Year.i=1998)
+    EddyProc.C$sPlotFingerprintY('NEE_f', Year.i=1998)
+    EddyProc.C$sPlotDailySumsY('NEE_f','NEE_fsd', Year.i=1998)
     EddyProc.C$sPlotDiurnalCycleM('NEE_f', Month.i=1)
 
     #+++ Export gap filled data to standard data frame
@@ -240,20 +242,19 @@ attr(sEddyProc.example,'ex') <- function(){
     #+++ Quality flag vector (e.g. from applying ustar filter)
     EddyDataWithPosix.F <- cbind(EddyDataWithPosix.F, QF=rep(c(1,0,1,0,1,0,0,0,0,0),nrow(EddyData.F)/10))
     #+++ Step function vector to simulate e.g. high/low water table
-    EddyDataWithPosix.F <- cbind(EddyDataWithPosix.F, Step= ifelse(EddyData.F$Day < 200 | EddyData.F$Day > 250, 0, 1))
+    EddyDataWithPosix.F <- cbind(EddyDataWithPosix.F, Step=ifelse(EddyData.F$DoY < 200 | EddyData.F$DoY > 250, 0, 1))
 
     #+++ Initialize eddy processing class with more columns
-    EddyTest.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE', 'LE', 'H', 'Rg', 'Tair', 'Tsoil', 'Rh', 'VPD', 'QF', 'Step'))
+    EddyTest.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE', 'LE', 'H', 'Rg', 'Tair', 'Tsoil', 'rH', 'VPD', 'QF', 'Step'))
 
-    #+++ Gap fill variable after applying quality flag QF with (non-default) variables and limits
+    #+++ Gap fill variable with (non-default) variables and limits including preselection with quality flag QF 
     EddyTest.C$sMDSGapFill('LE', QFVar.s='QF', QFValue.n=0, V1.s='Rg', T1.n=30, V2.s='Tsoil', T2.n=2, 'Step', 0.1, Verbose.b=T)
 
-    #+++ Try the gap filling subroutines individually for different window sizes and up to five variables and limits
-    EddyTest.C$sFillInit('NEE') #Initalize to fill variable 'NEE'
-    Result_Step1.F <- EddyTest.C$sFillLUT(3, 'Rg',50, 'Rh',15, 'Tair',2.5, 'Tsoil',1, 'Step',0.1)
+    #+++ Use individual gap filling subroutines with different window sizes and up to five variables and limits
+    EddyTest.C$sFillInit('NEE') #Initalize 'NEE' as variable to fill
+    Result_Step1.F <- EddyTest.C$sFillLUT(3, 'Rg',50, 'rH',30, 'Tair',2.5, 'Tsoil',2, 'Step',0.5)
     Result_Step2.F <- EddyTest.C$sFillLUT(6, 'Tair',2.5, 'VPD',3, 'Step',0.1)
     Result_Step3.F <- EddyTest.C$sFillMDC(3)
-    EddyTest.C$sPlotHHFluxes('VAR_f') #Individual fill result columns are called 'VAR_...'
-
+    EddyTest.C$sPlotHHFluxesY('VAR_fall', Year.i=1998) #Individual fill result columns are called 'VAR_...'
   }
 }

R/EddyGapfilling.R

@@ -51,13 +51,13 @@ sEddyProc$methods(
     ## VAR\emph{_fqc} - Quality flag assigned depending on gap filling method and window length (1 = most reliable, 2 = medium, 3 = least reliable)
     lTEMP <- data.frame(
       VAR_orig=Var.V.n               # Original values of variable VAR used for gap filling
-      ,VAR_f=NA                      # Original values and filled gaps
-      ,VAR_fall=NA                   # All values considered as gaps (for uncertainty estimates)
-      ,VAR_fnum=NA                   # Number of datapoints used for gap-filling
-      ,VAR_fsd=NA                    # Standard deviation of data points used for filling
-      ,VAR_fmeth=NA                  # Method used for gap filling
-      ,VAR_fwin= NA                  # Full window length used for gap filling
-      ,VAR_fqc= NA                   # Quality flag assigned depending on gap filling method and window length
+      ,VAR_f=NA_real_                # Original values and filled gaps
+      ,VAR_fall=NA_real_             # All values considered as gaps (for uncertainty estimates)
+      ,VAR_fnum=NA_real_             # Number of datapoints used for gap-filling
+      ,VAR_fsd=NA_real_              # Standard deviation of data points used for filling
+      ,VAR_fmeth=NA_real_            # Method used for gap filling
+      ,VAR_fwin=NA_real_             # Full window length used for gap filling
+      ,VAR_fqc=NA_real_              # Quality flag assigned depending on gap filling method and window length
     )
     sTEMP <<- data.frame(c(sTEMP, lTEMP))
 
@@ -99,7 +99,7 @@ sEddyProc$methods(
     #! Attention: For performance reasons, gap filled values and properties are first written to single variables and local matrix lGF.M
     #! (rather than changing single values in sTEMP which copies the data frame each time!)
     #! Improved algorithm speed by more than a factor of 10 (maybe even 100...)
-    lGF.M <- matrix(NA, nrow=0, ncol=7, dimnames=list(NULL,c('index','mean','fnum','fsd','fmeth','fwin','fqc')))
+    lGF.M <- matrix(NA_real_, nrow=0, ncol=7, dimnames=list(NULL,c('index','mean','fnum','fsd','fmeth','fwin','fqc')))
 
     # Check if sTEMP has been initialized with new VAR_ columns
     if( !exists('VAR_f', sTEMP) )
@@ -167,7 +167,7 @@ sEddyProc$methods(
 
           #Set window size and quality flag
           lVAR_fwin.n  <- 2*WinDays.i                      #! Full window length, congruent with MR PV-Wave, in paper single window sizes stated
-          lVAR_fmeth.n <- NA; lVAR_fqc.n <- NA;
+          lVAR_fmeth.n <- NA_real_; lVAR_fqc.n <- NA_real_;
           if( V1.s != 'none' && V2.s != 'none' && V3.s != 'none') { #Three conditions
             lVAR_fmeth.n <- 1
             if( lVAR_fwin.n <= 14 ) lVAR_fqc.n <- 1        #! Limit '14' congruent with MR PV-Wave, in paper different limit of '28' (stated as single window size of 14 days)
@@ -217,7 +217,7 @@ sEddyProc$methods(
     #! Attention: For performance reasons, gap filled values and properties are first written to single variables and local matrix lGF.M
     #! (rather than changing single values in sTEMP which copies the data frame each time!)
     #! Improved algorithm speed by more than a factor of 10 (maybe even 100...)
-    lGF.M <- matrix(NA, nrow=0, ncol=7, dimnames=list(NULL,c('index','mean','fnum','fsd','fmeth','fwin','fqc')))
+    lGF.M <- matrix(NA_real_, nrow=0, ncol=7, dimnames=list(NULL,c('index','mean','fnum','fsd','fmeth','fwin','fqc')))
 
     # Determine gap positions
     ToBeFilled.V.i <- which(is.na(sTEMP$VAR_fall))

R/FileHandling.R

@@ -60,29 +60,38 @@ attr(fLoadTXTIntoDataframe, 'ex') <- function() {
 
 fLoadFluxNCIntoDataframe <- function(
   ##title<<
-  ## Load specified variables and time stamp information from NetCDF file
+  ## Load NetCDF file
   ##description<<
-  ## The time stamp information needs to be provided as variables 'year', 'month', 'day', 'hour' (Fluxnet BGI format).
+  ## Load specified variables and time stamp information from NetCDF file in Fluxnet BGI format.
+  ## The time stamp information needs to be provided as variables 'year', 'month', 'day', 'hour'.
   VarList.V.s           ##<< Vector of variables to be read in 
   ,FileName.s           ##<< File name             
-  ,Dir.s=''             ##<< Directory      
+  ,Dir.s=''             ##<< Directory
+  ,NcPackage.s='ncdf4'  ##<< Name of R NetCDF package (implemented for 'RNetCDF' and 'ncdf4')
   ) 
   ##author<<
   ## AMM, KS
   # TEST: FileName.s <- 'Example_DE-Tha.1996.1998.hourly.nc'; Dir.s <- 'inst/MDSdata'; VarList.V.s <- c('NEE', 'Rg', 'rH', 'Tair', 'NEE_f')
+  # TEST: NcPackage.s <- 'ncdf4'
 {
+  # Check for R NetCDF packages
+  if( NcPackage.s=='ncdf4' ) { suppressWarnings(require(ncdf4))  
+  } else if( NcPackage.s=='RNetCDF' ) { suppressWarnings(require(RNetCDF))  
+  } else {
+    stop(CallFunction.s, ':::fLoadFluxNCIntoDataframe::: Required package \'RNetCDF\' or \'ncdf4\' could not be loaded!') }
+
   # Read in time variables
-  Data.F <- fAddNCFVar(NULL, 'year', FileName.s, Dir.s, 'fLoadFluxNCIntoDataframe')
-  Data.F <- fAddNCFVar(Data.F, 'month', FileName.s, Dir.s, 'fLoadFluxNCIntoDataframe')
-  Data.F <- fAddNCFVar(Data.F, 'day', FileName.s, Dir.s, 'fLoadFluxNCIntoDataframe')
-  Data.F <- fAddNCFVar(Data.F, 'hour', FileName.s, Dir.s, 'fLoadFluxNCIntoDataframe')
+  Data.F <- fAddNCFVar(NULL, 'year', FileName.s, Dir.s, NcPackage.s, 'fLoadFluxNCIntoDataframe')
+  Data.F <- fAddNCFVar(Data.F, 'month', FileName.s, Dir.s, NcPackage.s, 'fLoadFluxNCIntoDataframe')
+  Data.F <- fAddNCFVar(Data.F, 'day', FileName.s, Dir.s, NcPackage.s, 'fLoadFluxNCIntoDataframe')
+  Data.F <- fAddNCFVar(Data.F, 'hour', FileName.s, Dir.s, NcPackage.s, 'fLoadFluxNCIntoDataframe')
 
   # Convert time format to POSIX
   Data.F <- fConvertTimeToPosix(Data.F, 'YMDH', Year.s = 'year', Month.s='month', Day.s = 'day', Hour.s = 'hour')
 
   # Read in variables from a given list of needed variables 
   for (i in 1: length(VarList.V.s))  {
-    Data.F <- fAddNCFVar(Data.F, VarList.V.s[i], FileName.s, Dir.s, 'fLoadFluxNCIntoDataframe')
+    Data.F <- fAddNCFVar(Data.F, VarList.V.s[i], FileName.s, Dir.s, NcPackage.s, 'fLoadFluxNCIntoDataframe')
   }
   message('Loaded BGI Fluxnet NC file: ', FileName.s, ' with the following headers:')
   message('*** ', paste(colnames(Data.F), '(', as.character(lapply(Data.F, attr, which='units')), ')', collapse=' ', sep=''))
@@ -106,21 +115,17 @@ fAddNCFVar <- function(
   ,Var.s                ##<< Variable name
   ,FileName.s           ##<< NetCDF file name
   ,Dir.s                ##<< Directory
+  ,NcPackage.s          ##<< Name of R NetCDF package
   ,CallFunction.s=''    ##<< Name of function called from
   )
   ##author<<
   ## AMM, KS
   # TEST: Data.F <- NULL; Var.s <- 'NEE'; FileName.s <- 'Example_DE-Tha.1996.1998.hourly.nc'; Dir.s <- 'inst/MDSdata'
+  # TEST: NcPackage.s <- 'ncdf4'
 {
   InputNCF.s <- fSetFile(FileName.s, Dir.s, T, 'fAddNCFVar')
 
-  RNetCDF.b <- suppressWarnings(require(RNetCDF))
-  ncdf.b    <- suppressWarnings(require(ncdf4))
-
-  if ( !RNetCDF.b && !ncdf.b )
-      stop(CallFunction.s, ':::fAddNCFVar::: Required package RNetCDF or ncdf could not be loaded!') # for handling BGI Fluxnet netcdf files
-
-  if( RNetCDF.b ) {
+  if( NcPackage.s=='RNetCDF' ) {
     NCFile.C <- open.nc(InputNCF.s)
     tryCatch({
       NewCol.F <- data.frame(var.get.nc(NCFile.C, Var.s))
@@ -134,18 +139,14 @@ fAddNCFVar <- function(
     }, 
              finally = close.nc(NCFile.C)
     )
-  } else if( ncdf.b ) {
-#    stop('not implemented')
+  } else if( NcPackage.s=='ncdf4' ) {
     NCFile.C <- nc_open(InputNCF.s, write=FALSE, readunlim=TRUE, verbose=FALSE)
-      tryCatch({
-        NewCol.F <- data.frame(ncvar_get(NCFile.C, Var.s))
-        names(NewCol.F)[[1]] <- Var.s
-        attr(NewCol.F[[1]], 'varnames') <- Var.s
-        attr(NewCol.F[[1]], 'units') <- ncatt_get(NCFile.C, Var.s, 'units')
-
-      # Use c() instead of cbind() to be able to bind dataframe Data.F even if empty
+    tryCatch({
+      NewCol.F <- data.frame(ncvar_get(NCFile.C, Var.s))
+      names(NewCol.F)[[1]] <- Var.s
+      attr(NewCol.F[[1]], 'varnames') <- Var.s
+      attr(NewCol.F[[1]], 'units') <- ncatt_get(NCFile.C, Var.s, 'units')$value
       Data.F <- data.frame(c(Data.F, NewCol.F))
-      #attr(Data.F[[1]], 'units')
     }, 
              finally = nc_close(NCFile.C)
     )
@@ -164,88 +165,44 @@ fAddNCFVar <- function(
 
 fWriteDataframeToFile <- function(
   ##title<<
-  ## Write data frame to ASCII or NetCDF file
+  ## Write data frame to ASCII tab-separated text file
   Data.F                ##<< Data frame     
   ,FileName.s           ##<< File base name
   ,Dir.s=''             ##<< Directory
-  ,FileType.s='txt'     ##<< File output type
 )
   ##author<<
   ## AMM, KS
   ##details<<
-  ## 'txt' - for tab delimited text file with header and unit row
-  ## 'nc' - for very simple NetCDF file of numeric columns
-  ## With missing values flagged as -9999.0
-  # !!!TODO: NC file output: add non-numeric columns and units
-  # TEST: Data.F <- EddyData.F; BaseName.s <- 'OutputTest'; FileType.s='nc'; Dir.s <- 'data'; OutputName.s='none';
+  ## Missing values are flagged as -9999.0
+  # TEST: Data.F <- EddyData.F; FileName.s='none'; Dir.s <- 'data'; 
 {
   # Set file name
   OutputFile.s <- fSetFile(FileName.s, Dir.s, F, 'fWriteDataframeToFile')
 
   # Convert NAs to gap flag
   Data.F <- fConvertNAsToGap(Data.F)
 
-  # Write data to files
-  if( FileType.s=='txt') {
-    # Write tab delimited file
-    # supressWarnings()
-    Lines.V.s <- vector(mode='character', length = 2)
-    Lines.V.s[1] <- paste(colnames(Data.F), collapse='\t')
-    Lines.V.s[2] <- paste(as.character(lapply(Data.F, attr, which='units')), collapse='\t')
-    Lines.V.s[2] <- gsub('NULL', '--', Lines.V.s[2])
-    write(Lines.V.s, file=OutputFile.s, append=F)
-    write.table(format(Data.F, digits=5, drop0trailing=T, trim=T), file=OutputFile.s, col.names=F, row.names=F, sep='\t', quote=F, append=T)
-    message('Wrote tab separated textfile: ', OutputFile.s)
-
-  } else if( FileType.s=='nc') {
-    # Write NetCDF file
-    RNetCDF.b <- suppressWarnings(require(RNetCDF))
-    ncdf.b    <- suppressWarnings(require(ncdf4))
-
-    if ( !RNetCDF.b && !ncdf.b )
-      stop(CallFunction.s, ':::fWriteDataframeToFile::: Required package RNetCDF or ncdf could not be loaded!') # for handling BGI Fluxnet netcdf files
-
-    if( RNetCDF.b ) {
-      NCFile.C <- create.nc(OutputFile.s, clobber=T, large=T, prefill=F)
-      tryCatch({
-        dim.def.nc(NCFile.C, 'time', unlim=TRUE)
-        for (Var.i in 1:ncol(Data.F))  {
-          if( is.numeric(Data.F[,Var.i]) ) 
-          {
-            VarType.s <- 'NC_DOUBLE'
-            var.def.nc(NCFile.C, varname=names(Data.F)[Var.i], vartype=VarType.s, dimensions='time')
-            var.put.nc(NCFile.C, variable=names(Data.F)[Var.i], data=Data.F[,Var.i])
-            att.put.nc(NCFile.C, variable=names(Data.F)[Var.i], name='miss_val', type=VarType.s, value=-9999.0)
-          }
-          else next; #! Skips non-numeric columsn for now 
-        }
-      }, 
-               finally = close.nc(NCFile.C)
-      )
-    } else if( ncdf.b ) {
-      stop('!!! Error: not yet implemented !!!')
-      tryCatch({
-        NULL
-      },
-               finally = NULL
-      )
-    } else {
-      stop(CallFunction.s, ':::fWriteDataframeToFile::: NC files could not be opened!')
-    }
-    message('Wrote numeric columns to nc file: ', OutputFile.s)
-  }
+  # Write tab delimited file
+  Lines.V.s <- vector(mode='character', length = 2)
+  Lines.V.s[1] <- paste(colnames(Data.F), collapse='\t')
+  Lines.V.s[1] <- gsub('DateTime', 'Date Time', Lines.V.s[1]) #POSIX column
+  Lines.V.s[2] <- paste(as.character(lapply(Data.F, attr, which='units')), collapse='\t')
+  Lines.V.s[2] <- gsub('NULL', '-', Lines.V.s[2])
+  Lines.V.s[2] <- gsub('DateTime', 'Date Time', Lines.V.s[2])  #POSIX column
+  write(Lines.V.s, file=OutputFile.s, append=F)
+  write.table(format(Data.F, digits=5, drop0trailing=T, trim=T), file=OutputFile.s, col.names=F, row.names=F, sep='\t', quote=F, append=T)
+  message('Wrote tab separated textfile: ', OutputFile.s)
 
-##value<<
+  ##value<<
   ## Output of data frame written to file of specified type.
 }
 
 attr(fWriteDataframeToFile, 'ex') <- function() {
   # Example code
+  if (FALSE) { #Example code, do not always execute (e.g. on package installation)
   if( file.exists('data/Example_DETha98.txt') ) {
-    EddyData.F <- fLoadTXTIntoDataframe('Example_DETha98.txt','data')
-    if (FALSE) { #Example code, do not always execute (e.g. on package installation)
-      fWriteDataframeToFile(EddyData.F, 'OutputTest', 'out')
-      fWriteDataframeToFile(EddyData.F, 'OutputTest', 'out', 'nc')
+      EddyData.F <- fLoadTXTIntoDataframe('Example_DETha98.txt','data')
+      fWriteDataframeToFile(EddyData.F, 'OutputTest.txt', 'out')
     }
   }
 }
@@ -309,29 +266,30 @@ fSetFile <- function(
 )
   ##author<<
   ## AMM
-  # TEST: Dir.s <- 'data'; FileName.s <- 'test'; FileName.s <- 'Example_DETha98.txt'; IO.b <- T
+  # TEST: Dir.s <- 'data'; FileName.s <- 'Example_DETha98.txt'; IO.b <- T; CallFunction.s <- 'test'
 {
   # Check if string for directory provided
   Dir.b <- fCheckValString(Dir.s)
 
   # Check if directory exists
-  if ( Dir.b && !file.exists(Dir.s) && IO.b )
-    stop(CallFunction.s, '::: Directory does not exist: ', Dir.s)
+  if ( Dir.b && (file.access(Dir.s, mode=4) != 0) && IO.b )
+    stop(CallFunction.s, ':::fSetFile::: Directory does not exist: ', Dir.s)
 
   # Make directory if mode is output
-  if( Dir.b && !file.exists(Dir.s) && !IO.b ) {
+  if( Dir.b && (file.access(Dir.s, mode=0) != 0) && !IO.b ) {
     dir.create(Dir.s)
-    if( !file.exists(Dir.s) )
-      stop(CallFunction.s, '::: Directory could not be created: ', Dir.s)
+    message(CallFunction.s, ':::fSetFile::: Directory created: ', Dir.s)
+    if( file.access(Dir.s, mode=2) != 0 )
+      stop(CallFunction.s, ':::fSetFile::: Directory could not be created: ', Dir.s)
   }
 
   # Set file name accordingly
   File.s <- if( Dir.b ) {  paste(Dir.s, '/', FileName.s, sep='')
   } else { FileName.s }
 
   # If input file, check if file exists
-  if ( IO.b && !file.exists(File.s))
-    stop(CallFunction.s, '::: File does not exist: ', File.s)
+  if ( IO.b && (file.access(File.s, mode=4) != 0) )
+    stop(CallFunction.s, ':::fSetFile::: File does not exist: ', File.s)
 
   File.s
   ##value<<