#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#+++ Developers' R script +++
#+++ Various test routines to develop and test code and processing for various datasets
#+++ --> run from package directory
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Author: AMM
Develop.b <- F #True if in development mode: Load individual scripts, if false test REddyProc as package
ShortTest.b <- F #Short test only
LongTest.b <- F #True if in intensive test mode including NC files, all plots and ...
#LongTest.b <- T #True if in intensive test mode including NC files, all plots and ...
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Load package or source files (development mode)
if (Develop.b) {
# Source settings for R environment and standard functions
source('inst/develop/setREnvir.R')
# Source file and data handling scripts
source('R/DataFunctions.R')
source('R/FileHandling.R')
# Source geo functions
source('R/GeoFunctions.R')
# Initialization of sEddyProc class
source('R/Eddy.R')
# sEddyProc methods
source('R/EddyGapfilling.R')
source('R/EddyPlotting.R')
source('R/EddyPartitioning.R')
# sEddyProc method still in development
source('inst/develop/EddyFiltering.R')
} else {
# Source settings for R environment and standard functions
source('inst/develop/setREnvir.R')
# Generate package from local path
### system('R CMD INSTALL --build --html --library = /Library/Frameworks/R.framework/Versions/current/Resources/library ../REddyProc')
### system('R CMD INSTALL --build --html ../REddyProc')
# Test installation of package from various sources
# Requires restart of R console afterwards!
### R-Forge: install.packages("REddyProc", repos = "http://R-Forge.R-project.org", type = "source")
### R-Forge after manual download: install.packages('../_FromRForge/REddyProc_0.5.tar.gz', repos = NULL)
### Self generated package: install.packages('REddyProc_0.5.tgz', repos = NULL)
require('REddyProc')
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Load example data directly from package or (if not available) from txt file
#data('Example_DETha98', package = 'REddyProc') # superseedeby by LazyData: true
EddyData.F <- if (exists("Example_DETha98") ) {
# example is available only when the package is loaded
Example_DETha98
} else {
# if the package is not loaded, assume to be some directory relative to the example
if (file.exists('../examples/Example_DETha98.txt') ) {
Example_DETha98 <- suppressMessages(fLoadTXTIntoDataframe('Example_DETha98.txt','../examples'))
} else {
message('Unit test directory: ', getwd())
message('Workspace: ', ls())
stop('test_sEddyProc.R::: Example data could not be loaded.')
}
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Run short test
if (ShortTest.b ) {
EddyDataWithPosix.F <- fConvertTimeToPosix(EddyData.F, 'YDH', Year = 'Year', Day = 'DoY', Hour = 'Hour')
EPTha.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE', 'Rg', 'Tair', 'VPD', 'Ustar'))
EPTha.C$sMDSGapFill('NEE')
EPTha.C$sMDSGapFill('Tair',FillAll = F)
EPTha.C$setLocationInfo(LatDeg = 51.0, LongDeg = 13.6, TimeZoneHour = 1.0)
EPTha.C$sMRFluxPartition()
View(EPTha.C$sTEMP)
stop('No error but short test only.')
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Run unit tests
require(testthat) #Packages to source for unit tests
if (Develop.b) {
test_dir('inst/tests') # works only if R scripts are sourced
} else {
test_package('REddyProc') # works only if package is installed (not necessarily loaded)
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Testing of ustar filtering
if (FALSE ) { #Short example from above with ustar filtering
EddyDataWithPosix.F <- fConvertTimeToPosix(EddyData.F, 'YDH', Year.s = 'Year', Day.s = 'DoY', Hour.s = 'Hour')
EPTha.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE', 'Rg', 'Tair', 'VPD', 'Ustar'))
EPTha.C$sUstarMM() # Apply filter
EPTha.C$sMDSGapFill('NEE','UstarMM_fqc',0) #Use ustar flag
EPTha.C$sMDSGapFill('Tair',FillAll = F)
EPTha.C$setLocationInfo(LatDeg = 51.0, LongDeg = 13.6, TimeZoneHour = 1.0)
EPTha.C$sMRFluxPartition()
View(EPTha.C$sTEMP)
stop('No error but short test only.')
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Data handling of standard example data
# Add dummy quality flag for tests
EddyData.F <- cbind(EddyData.F, QF = structure(rep(c(1,0,1,0,1,0,0,0,0,0),nrow(EddyData.F)/10), units = 'dummy_unit'))
# Test calculation of VPD
EddyData.F$VPDnew <- fCalcVPDfromRHandTair(EddyData.F$rH, EddyData.F$Tair)
# Add POSIX time stamp
EddyDataWithPosix.F <- fConvertTimeToPosix(EddyData.F, 'YDH', Year = 'Year', Day = 'DoY', Hour = 'Hour')
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Load netcdf datasets in BGI format
if (LongTest.b) {
# Load NC file with multiple years (upload includes time conversion)
lVar.V.s <- c('NEE', 'Rg', 'Tair', 'VPD', 'NEE_f', 'NEE_fmet', 'NEE_fwin', 'NEE_fn', 'NEE_fs', 'NEE_fqc', 'NEE_fqcOK') #!!! Attention NEE_fqcOK or NEE_fqcok, both exists
#EddyNCData.F <- fLoadFluxNCIntoDataframe(lVar.V.s, getExamplePath('Example_DE-Tha.1996.1998.hourly_selVars.nc',TRUE))
}
# Run loop over all (site) files in BGI Fluxnet data directory
if (T == F) {
SiteFile.V.s <- fInitFilesDir(DirFluxnet.s, 'hourly.nc')
SiteName.V.s <- fStripFileExtension(SiteFile.V.s)
for (Site.i in 1:length(SiteName.V.s)) {
message(paste('Handling site file ', Site.i, ': \'', SiteName.V.s[Site.i],'\'', sep = ''))
#...
}
#EddyBGINCData.F <- fLoadFluxNCIntoDataframe(lVar.V.s, 'DE-Tha.1996.2006.hourly.nc', DirFluxnet.s)
# EddyBGINCData.F <- fLoadFluxNCIntoDataframe(lVar.V.s, 'DE-Tha.1996.2006.hourly.nc', DirFluxnet.s,'RNetCDF')
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Write data to files
if (LongTest.b) {
fWriteDataframeToFile(EddyDataWithPosix.F, 'DE-Tha-Data.txt', 'out')
}
#Produce new ascii test files from BGI netcdf fluxnet files
if (FALSE) {
lVar2.V.s <- c('NEE', 'LE', 'H', 'Rg', 'VPD', 'rH', 'Tair', 'Tsoil_f', 'julday')
#Example.F <- fLoadFluxNCIntoDataframe(lVar2.V.s, getExamplePath('Example_DE-Tha.1996.1998.hourly_selVars.nc',TRUE))
#Example.F <- fLoadFluxNCIntoDataframe(lVar2.V.s, getExamplePath('Example_DE-Tha.1996.1998.hourly_selVars.nc'), 'RNetCDF')
Example.F$Year <- as.numeric(format(Example.F$DateTime, '%Y'))
Example.F$Month <- as.numeric(format(Example.F$DateTime, '%m'))
Example.F$DoY <- as.numeric(format(Example.F$DateTime, '%j'))
Example.F$Hour <- as.numeric(format(Example.F$DateTime, '%H')) + as.numeric(format(Example.F$DateTime, '%M'))/60
colnames(Example.F)[colnames(Example.F) == 'Tsoil_f'] <- 'Tsoil'
fWriteDataframeToFile(Example.F, 'DE-Tha.1996.1998.txt','out')
# Try to reload data file
Eddy3Years.F <- fLoadTXTIntoDataframe('DE-Tha.1996.1998.txt','out')
Eddy3Years.F <- fConvertTimeToPosix(Eddy3Years.F, 'YDH', Year.s = 'Year', Day.s = 'DoY', Hour.s = 'Hour')
#fCheckHHTimeSeries(Eddy3Years.F$DateTime, DTS = 48)
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Load data into REddyProc class instances
# Standard dataset
EPTha.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE', 'QF', 'Rg', 'Tair', 'VPD'))
# Dataset with pseudo hourly data
EPThaH.C <- sEddyProc$new('DE-Tha-H', EddyDataWithPosix.F[c(F,T),], c('NEE', 'QF', 'Rg', 'Tair', 'VPD'), DTS.n = 24)
# Subset of dataset
EPThaS.C <- sEddyProc$new('DE-Tha-S', EddyDataWithPosix.F[(4321:8640),], c('NEE', 'QF', 'Rg', 'Tair', 'VPD'))
# Limited dataset with NEE only
EPThaL1.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE'))
EPThaL2.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE', 'Rg'))
EPThaL3.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE', 'Rg', 'Tair'))
# Work with (subsets of) multiple years
if (LongTest.b) {
# All three years of data
EPThaNC.C <- sEddyProc$new('DE-Tha', EddyNCData.F, lVar.V.s)
# Single year 1996
EPThaNC96.C <- sEddyProc$new('DE-Tha', EddyNCData.F[1:17568,], lVar.V.s)
# Subset of half 1996 to half of 1997
EPThaNCsub.C <- sEddyProc$new('DE-Tha', EddyNCData.F[8737:(17568+9600),], lVar.V.s)
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Fill gaps with MDS algorithm
#system.time(...)
EPTha.C$sMDSGapFill('NEE','QF', 0, FillAll = T, isVerbose = T)
EPTha.C$sMDSGapFill('NEE', V1.s = 'none', FillAll = T, isVerbose = T)
EPThaH.C$sMDSGapFill('NEE','QF', 0, FillAll = T, isVerbose = T)
EPThaS.C$sMDSGapFill('NEE','QF', 0, FillAll = T, isVerbose = T)
EPThaL1.C$sMDSGapFill('NEE', V1.s = 'none', T1.n = NA_real_, V2.s = 'none', T2.n = NA_real_, V3.s = 'none', T3.n = NA_real_, FillAll = T, isVerbose = T)
EPThaL1.C$sMDSGapFill('NEE', V1.s = 'none', FillAll = T, isVerbose = T)
EPThaL1.C$sMDSGapFill('NEE', FillAll = T, isVerbose = T)
EPThaL2.C$sMDSGapFill('NEE', FillAll = T, isVerbose = T)
EPThaL3.C$sMDSGapFill('NEE', FillAll = T, isVerbose = T)
# Fill also other variables
EPTha.C$sMDSGapFill('Rg', FillAll = T, isVerbose = T)
EPTha.C$sMDSGapFill('VPD', FillAll = T, isVerbose = T)
EPTha.C$sMDSGapFill('Tair', FillAll = T, isVerbose = T)
# Fill multiple years
if (LongTest.b) {
EPThaNC.C$sMDSGapFill('NEE_f', 'NEE_fqc', 0, FillAll = T, isVerbose = T)
EPThaNC96.C$sMDSGapFill('NEE_f', 'NEE_fqc', 0, FillAll = T, isVerbose = T)
EPThaNCsub.C$sMDSGapFill('NEE_f', 'NEE_fqc', 0, FillAll = T, isVerbose = T)
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Export processing results
ThaFilled.F <- EPTha.C$sExportResults()
fWriteDataframeToFile(cbind(EddyDataWithPosix.F,ThaFilled.F), 'DE-Tha-ResultsTest.txt', 'out')
if (LongTest.b) {
ThaHFilled.F <- EPThaH.C$sExportResults()
ThaSFilled.F <- EPThaS.C$sExportResults()
ThaNCFilled.F <- EPThaNC.C$sExportResults()
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Test documention, e.g. information on reference class methods
# sEddyProc$help(sMDSGapFill)
#Show function code
# fix('sEddyProc.example')
# fix('fSetQF')
# fix('EPTha.C')
# fix('sEddyProc$sMDSGapFill') or fix('EPTha.C$sMDSGapFill') --> not possible to access S4 functions... :-(
# Show data in reference class
# EPTha.C$sPrintData()
# View(EPTha.C$sDATA)
# View(EPTha.C$sTEMP)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++]
# Plotting of single years/month
fPlots <- function(ClassName.s, Var.s, QFVar.s = 'none', QFValue.n = NA, VarUnc.s = 'none', Month.i = NA, Year.i = NA) {
if (is.na(Year.i)) { #All years/months and save to file
eval(parse(text = paste(ClassName.s, '$sPlotFingerprint(\'', Var.s, '\', QFVar.s = \'', QFVar.s, '\', QFValue.n = ', QFValue.n, ')', sep = '')))
eval(parse(text = paste(ClassName.s, '$sPlotHHFluxes(\'', Var.s, '\', QFVar.s = \'', QFVar.s, '\', QFValue.n = ', QFValue.n, ')', sep = '')))
eval(parse(text = paste(ClassName.s, '$sPlotDailySums(\'', Var.s, '\', VarUnc.s = \'', VarUnc.s, '\')', sep = '')))
eval(parse(text = paste(ClassName.s, '$sPlotDiurnalCycle(\'', Var.s, '\', QFVar.s = \'', QFVar.s, '\', QFValue.n = ', QFValue.n, ')', sep = '')))
} else { #Individual years/months
eval(parse(text = paste(ClassName.s, '$sPlotFingerprintY(\'', Var.s, '\', QFVar.s = \'',
QFVar.s, '\', QFValue.n = ', QFValue.n, ', Year.i = ', Year.i, ')', sep = '')))
eval(parse(text = paste(ClassName.s, '$sPlotHHFluxesY(\'', Var.s, '\', QFVar.s = \'',
QFVar.s, '\', QFValue.n = ', QFValue.n, ', Year.i = ', Year.i, ')', sep = '')))
eval(parse(text = paste(ClassName.s, '$sPlotDailySumsY(\'', Var.s, '\', VarUnc.s = \'',
VarUnc.s, '\', Year.i = ', Year.i, ')', sep = '')))
eval(parse(text = paste(ClassName.s, '$.sPlotDiurnalCycleM(\'', Var.s, '\', QFVar.s = \'',
QFVar.s, '\', QFValue.n = ', QFValue.n, ', Month.i = ', Month.i, ')', sep = '')))
}
}
if (LongTest.b) {
#Before filling
fPlots('EPTha.C', 'NEE', Month.i = 10, Year.i = 1998)
fPlots('EPThaH.C', 'NEE', Month.i = 10, Year.i = 1998)
fPlots('EPThaS.C', 'NEE', Month.i = 10, Year.i = 1998)
fPlots('EPThaNCsub.C', 'NEE', Month.i = 6, Year.i = 1996)
fPlots('EPThaNCsub.C', 'NEE', Month.i = 6, Year.i = 1997)
#After filling
fPlots('EPTha.C', 'NEE_f','NEE_fqc', 1, 'NEE_fsd', Month.i = 10, Year.i = 1998)
fPlots('EPThaH.C', 'NEE_f','NEE_fqc', 1, 'NEE_fsd', Month.i = 10, Year.i = 1998)
fPlots('EPThaS.C', 'NEE_f','NEE_fqc', 1, 'NEE_fsd', Month.i = 10, Year.i = 1998)
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Plotting of standard example, unfilled and filled columns
fPlots('EPTha.C', 'NEE')
fPlots('EPTha.C', 'NEE_f')
fPlots('EPTha.C', 'NEE_f', 'NEE_fqc', 1)
# Images for multiple year data
if (LongTest.b) {
# Multiple year plots with leap year
fPlots('EPThaNC.C', 'NEE')
fPlots('EPThaNC.C', 'NEE_f')
fPlots('EPThaNC.C', 'NEE_f', 'NEE_fqc', 1)
fPlots('EPThaNC96.C', 'NEE')
fPlots('EPThaNCsub.C', 'NEE')
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Re-set class on previously filled data for plotting, e.g. to avoid rerunning MDS each time...
if (LongTest.b) {
EPTha.C <- sEddyProc$new('DE-Tha', cbind(EddyDataWithPosix.F,ThaFilled.F), c('NEE','Rg', 'Tair', 'VPD', 'QF', 'NEE_f', 'NEE_fqc', 'NEE_fmeth', 'NEE_fwin', 'NEE_fsd', 'NEE_fnum'))
EPTha.C$sPlotDailySums('NEE_f','NEE_fsd')
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Compare data of old MDS filling with new filling
if (LongTest.b) {
# Standard filling - without ustar filtering
EPTha98.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE', 'QF', 'Rg', 'Tair', 'VPD'))
EPTha98.C$sMDSGapFill('NEE', FillAll = T, isVerbose = T)
# Load MDS output data from old PV-Wave online tool - without ustar filtering
MDSData.F <- fLoadTXTIntoDataframe('Example_DETha98_PVWave_DataSetafterGapfill.txt','inst/examples')
MDSData.F <- fConvertTimeToPosix(MDSData.F, 'YMDHM', Year.s = 'Year', Month.s = 'Month', Day.s = 'Day', Hour.s = 'Hour', Min.s = 'Minute')
# Plot difference between old and new MDS
plot(EPTha98.C$sTEMP$NEE_f ~ MDSData.F$NEE_f, col = rgb(0.4,0.4,0.4,alpha = 0.2), pch = 20, cex = 0.3)
# Load MDS output data from LaThuile BGI fluxnet files with ustar filtering
# Rename 'NEE' to have different column names for old and new processing
# (though in different data frames: original data in sDATA and results in sTEMP)
EddyNCData.F <- cbind(EddyNCData.F, NEEnew = EddyNCData.F$NEE)
EPThaNC.C <- sEddyProc$new('DE-Tha', EddyNCData.F, c('NEEnew', 'Rg', 'Tair', 'VPD', 'NEE_f', 'NEE_fqc', 'NEE_fmet', 'NEE_fwin', 'NEE_fs', 'NEE_fn'))
# Fill gaps - with ustar through NEE_fqc = 0
EPThaNC.C$sMDSGapFill('NEEnew', 'NEE_fqc', 0, FillAll = T, isVerbose = T)
# Plot difference between old and new MDS
plot(EPThaNC.C$sTEMP$NEEnew_f ~ EPThaNC.C$sDATA$NEE_f, col = rgb(0.4,0.4,0.4,alpha = 0.2), pch = 20, cex = 0.3)
# Plot combination of the two BUT this is mixing with and without ustar filtering...
plot(EPTha98.C$sTEMP$NEE_f ~ EPThaNC.C$sTEMP$NEEnew_f[35089:52608], col = rgb(0.4,0.4,0.4,alpha = 0.2), pch = 20, cex = 0.3)
plot(MDSData.F$NEE_f ~ EPThaNC.C$sDATA$NEE_f[35089:52608], col = rgb(0.4,0.4,0.4,alpha = 0.2), pch = 20, cex = 0.3)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.