bgctw / reddyproc Goto Github PK
View Code? Open in Web Editor NEWProcessing data from micrometeorological Eddy-Covariance systems
Processing data from micrometeorological Eddy-Covariance systems
Currently, night-time based flux partitioning estimates temperature sensitivity E0 for a shifting time window before adopting a single aggregated value. There is a constraint on temperature range before a window is considered for fitting the respiration-temperature relationship. By default temperature range must be at least TempRange.n = 5°Celsius. The rational: with little variation in temperature, there will not be a good and maybe spurious relationship with respiration.
At some sites, the night-time temperature range over seven days rarely changes this much, leading to an error message in night-time based flux partitioning. Although E0 estimates become less reliable with less temperature variation, the user should be able to try relaxing this criterion.
The criterion can already specified at functions fRegrE0fromShortTerm
and also at sEddyProc_sMRFluxPartition
with argument parsE0Regression
. But the user only gets the
warning that respiration-temperature relationship could not be fit at enough periods.
I suggest telling the user in the warning how he can retry with a relaxed temperature range constraint.
Currently, REddyProc checks for columns to be numeric during class initialization. If columns are non-numeric the column names that are non-numeric are reported by an error message.
In order to help the user to find the violations, find the first non-numeric value and report it together with the row in the error message.
Currently the ticks at the fingerprint plots are show at the center of the month. They should better appear at the beginning of the month.
Currently, the error message in case of missing Ustar column is misleading, indicating missing temperature variable.
Check for all required columns and report missing columns in error message.
Currently, Setting "isVerbose = FALSE" only suppresses the print of "Start daytime flux partitioning ..." but not "Estimating temperature sensitivity from night time NEE , 1, 3, 5, 7, 9, ...". It should suppress all progress output.
Currently, first the full model (including the VPD effect) is fitted, and only if this is not successful reduced versions are fitted.
If there is no large VPD in the dataset, one could omit the full model and start with model version neglecting the VPD.
Keenan 2019 proposed a modified daytime partitioning that accounts for inhibition of leaf respiration during light. It modifies GLPartitioning by estimating night-time respiration by basal respiration estimated at night instead of estimated from daytime LRC curve.
Currently, if there is a year where no respiration~temperature relationship could be established in any time window, the daytime partitioning quits with an error.
Replace this behavior by issuing a warning and use the mean temperature sensitivity of the other years for the year in question. Use a reasonable upper standard deviation: a bit higher than the a robust maximum: quantile(sdE0, 0.9) * 1.2
Currently, the flux partitioning is repeatedly called for different uStar scenarios and correctly warns that columns are replaced, because the uStar-related columns are overridden.
However, in this case this is intended behaviour and the warning should be avoided.
Hi,
I am trying to do a workaround for the behaviour that is not allowing me to store all the parameters of sGLFluxPartition when run over multiple suffixes. If that is not too problematic due to the internal structure, it would be better if it can be corrected in the first place. If I use suffixes c('uStar', 'U05', 'U50', 'U95', 'UNone'), I will get some of the columns with suffixes after sExportResults
, while some not. Some of them are identical for different runs, some of them not (PotRad, FP_Temp should not get suffixes as they are identical, also in sMRFluxPartition
). Those variables without suffix will be overwritten, I will get only the results for 'UNone' and Warning message is shown:
Warning message:
In EC$sGLFluxPartition(TempVar.s = paste(FP_temp, "_f", sep = ""), :
replacing existing output columnsFP_VARnight, FP_VARday, NEW_FP_Temp, NEW_FP_VPD, FP_RRef_Night, FP_qc, FP_dRecPar, FP_errorcode, FP_GPP2000, FP_OPT_VPD, FP_OPT_NoVPD, FP_k, FP_beta, FP_alpha, FP_RRef, FP_E0, FP_k_sd, FP_beta_sd, FP_alpha_sd, FP_RRef_sd, FP_E0_sd, Reco_DT_uStar, GPP_DT_uStar, Reco_DT_uStar_SD, GPP_DT_uStar_SD
When I go through the suffix loop, resave the object and run sGLFluxPartition
for 'uStar' suffix again I can check that content of these columns is not identical and thus should have respective suffix:
[1] "Reco_DT_uStar_SD" "GPP_DT_uStar_SD" "FP_VARnight" "FP_VARday" "FP_RRef_Night"
[6] "FP_qc" "FP_dRecPar" "FP_errorcode" "FP_GPP2000" "FP_k"
[11] "FP_beta" "FP_alpha" "FP_RRef" "FP_E0" "FP_k_sd"
[16] "FP_beta_sd" "FP_alpha_sd" "FP_RRef_sd" "FP_E0_sd"
Please ignore "Reco_DT_uStar_SD", "GPP_DT_uStar_SD". These are probably result of some random number generation during SD reestimation for the second time. Also the warning message is complaining about rewriting of variables with 'uStar' suffix which is correct behaviour. You can notice that warning message is also mentioning collumns that are not finally exported.
I would suggest to include suffix support for the above vector of variable names and relax the expectations for the unexported ones. I assume it will not affect the correct working of the method. Maybe you could also include suffixes in the internal collumns and get rid of warning()
whatsoever - replacing original content with almost identical one should not cause issues.
Thanks!
Hi,
I just noticed in the vignette that the condition if (nEff == 0)
is not strong enough when computing sdNEE
due to division by (nEff - 1)
. In case that nEff == 1
computation will produce sdNEE == Inf
.
Example:
aggDay[174, ]
In my eddy convince data, there is no Ustar data. How can i generate the half-hour GPP by using the half-hour 'Year', 'Doy', 'NEE', 'Rg' and 'Tair'.
Hi,
REddyProc is a great package. I wonder whether this package would realize the function of footprint modeling in Tovi in the future?
Thanks
When using sEddyProc_sPlotNEEVersusUStarForSeason
following warning is displayed with suggested code correction:
Warning message:
`filter_()` was deprecated in dplyr 0.7.0.
Please use `filter()` instead.
See vignette('programming') for more help
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
Seems to be connected with changes in tidyverse packages. The warning was displayed already in REddyProc_1.3.0 and is still present in REddyProc_1.3.1.
In order to compare to Fluxnet15, provide the results in format of the Fluxnet 2015 release
https://fluxnet.fluxdata.org/data/fluxnet2015-dataset/fullset-data-product/
Hi all, my flux tower data has 15 minutes resolution, so I will work to modify my REddyProc copy. Suggestions to start? I´m new at R and learning as fast as possible. Thanks. Leila
There is be a small bug that it throws error when arg probs
is assigned only one numeric value. I hope you can replicate the problem with this command:
(uStarRes <- EddyProc.C$sEstUstarThresholdDistribution(
nSample = 5L, probs = 0.5,
ctrlUstarEst.l = usControlUstarEst(isUsingCPTSeveralT = TRUE)))
Error:
Error in `[<-`(`*tmp*`, iInvalid, , value = NA_real_) :
(subscript) logical subscript too long
Actually it works if I repeat the probability value:
(uStarRes <- EddyProc.C$sEstUstarThresholdDistribution(
nSample = 5L, probs = c(0.5, 0.5),
ctrlUstarEst.l = usControlUstarEst(isUsingCPTSeveralT = TRUE)))
Thank you for your feedback! It was really helpful.
currently each user has to access the names of sGetUstarScenarios()
. Moreover, this gives warnings, if no uStar thresholds have been set yet. Better provide a method support a simpler handling.
Hi all,
I am new in Eddy flux data processing.
I working on AsiaFLux data.
The dataset is already in half-hourly format, the data frame looks like this:
head(BKS_Eddy02)
Year DoY Hour NEE LE H Rg Tair Tsoil rH VPD Ustar
1 NA NA NA micro mol m-2 s-1 Wm-2 Wm-2 Wm-2 degrees C degreesC % kPa ms-1
2 2001 1 15 -99999 -99999 -99999 0 24.278 25.18 96 0.1188 -99999
3 2001 1 45 -99999 -99999 -99999 0 24.31 25.138 96.075 0.1169 -99999
4 2001 1 115 1.895 -99999 -99999 0 24.245 25.087 96.25 0.1112 -99999
5 2001 1 145 -99999 -99999 -99999 0 24.36 25 96.225 0.1127 -99999
6 2001 1 215 -25.448 -99999 -99999 0 25.035 24.97 94.475 0.1726 -99999
When I tried to convert the date and time into Posix using the EdyyProc package, I got an error message.
BKS_Eddy02$DateTime <- fConvertTimeToPosix(BKS_Eddy02, 'YDH', Year = 'Year', Day = 'DoY', Hour = 'Hour')
Error in if (sum(is_invalid_doy) > 0) { :
missing value where TRUE/FALSE needed
How should I treat my data? The time format seems different from the other data from FLuxnet or with the example data.
Thank you in advance
Currently, the user has to provide the association between seasons and uStar and has to provide suffixes for flux partitioning. It would be nice, to support setting uStar scenarios once and next invoke simple commands to process across all the scenarios.
The Lasslop daytime partitioning of GPP currently relies on smoothing parameters using a Gaussian Process regression using package mlegp. However, mlegp is being archived from CRAN threatening to also removing REddyProc from CRAN.
Hopefully mlegp authors will fix this issue. Otherwise
In the longer-time do the smoothing using some other method.
In the short-time we need to replace the strong dependency to mlegp by a weaker dependency (in Suggests) and require the user to install a maybe old version of mlegp by hand.
email from CRAN maintainers: "
We have asked for corrections to the errors showing at
https://www.stats.ox.ac.uk/pub/bdr/donttest/FCPS.out
https://www.stats.ox.ac.uk/pub/bdr/donttest/penalizedSVM.out
with no response from the mleqp and densityClust maintainers where the errors originate. That for mlegp even shous in Rdev CMD check --as-cran:
Thus those packages are now scheduled for archival on 2022-03-20, and archiving them will necessitate also archiving their strong reverse dependencies
REddyProc REddyProcNCDF TAG c060 monocle penalizedSVM
if they still require one of these packages by that date.
"
e.g. VPD computation is better maintained in the bigleaf package than in REddyProc
NAs for seasonal uStar Threshold estimates can be introduced during the bootstrap in sEstUstarThresholdDistribution if the number of valid bootstraps is lower than usControlUstarEst()$minValidBootProp.
When extracting seasonal values in usGetSeasonalSeasonUStarMappingFromDistributionResult, these NA values should be replaced by the annual estimates after a warning.
REddyProc sometimes fails or provides strange results, if unreasonable NEE uncertainty estimates are provided. E.g. at some sites, NEE is repeatedly at the numerically exact same value at similar daytime, leasing to a zero uncertainty estimates. This in turn leads to non-finite values of cost-functions, e.g. in LRC fitting, and further upstream problems, whose cause is hard very hard to track down.
Therefore, REddyProc should early on detect such problematic sdNEE records, warn the user and do something about it before further processing.
I suggest the following
The method will inform user about missing Tair_f or Rg_f but not about missing VPD_f.
Hi,
I was comparing with my results using base R functions and could not arrive to those presented for daily aggregation in vignette. After some testing I noticed that %>% filter(NEE_uStar_fqc == 0)
is missing in the pipe, i.e. all daily values of NEE_uStar_fsd
are used for aggregation without filtering. I confirmed when I removed the filtering in my code. Unfortunately, the solution is not so easy as including above filter in the pipe will also filter nEff
that will again affect the uncertainty results.
For historical development reasons, currently, the type of an argument is encoded in its name.
This is confusing to the user.
For the most important functions, omit the type suffix from argument names. But at the same time make sure that code using the old argument names still runs, and gives a warning of changed argument names.
Hello,
I tried several things but I still get NA uStar thresholds when using CPT. I am running the newest R version 3.5.1 (2018-07-02), the newest RStudio Version 1.1.456 and CRAN REddyProc 1.1.5. I tried to update the tidyverse packages as suggested on GitHub. I was not able to install REddyProc from GitHub through devtools::install_github("bgctw/REddyProc"). This seems to be unrelated issue that devtools currently do not recognize Rtools produced by Rtools35.exe. Here is a reproducible example with some console outputs:
library(REddyProc)
fileName <- getExamplePath('Example_DETha98.txt', isTryDownload = TRUE)
EddyData.F <- if (length(fileName)) fLoadTXTIntoDataframe(fileName) else Example_DETha98
EddyDataWithPosix.F <- fConvertTimeToPosix(EddyData.F, 'YDH',Year.s = 'Year',Day.s = 'DoY',Hour.s = 'Hour')
EddyProc.C <- sEddyProc$new('DE-Tha', EddyDataWithPosix.F, c('NEE','Rg','Tair','VPD', 'Ustar'))
EddyProc.C$sEstUstarThreshold(ctrlUstarEst.l = usControlUstarEst(isUsingCPTSeveralT = TRUE))$uStarTh
With output:
aggregationMode seasonYear season uStar
1 single NA NA
2 year 1998 NA
3 season 1998 1998001 NA
4 season 1998 1998003 NA
5 season 1998 1998006 NA
6 season 1998 1998009 NA
7 season 1998 1998012 NA
and alternative:
season_factor <- usCreateSeasonFactorMonthWithinYear( EddyDataWithPosix.F$DateTime - 15 * 60)
table(season_factor)
EddyProc.C$sEstUstarThreshold(ctrlUstarEst.l = usControlUstarEst( isUsingCPTSeveralT = TRUE), seasonFactor.v = season_factor)$uStarTh
with similar results:
aggregationMode seasonYear season uStar
1 single NA NA
2 year 1998 NA
3 season 1998 1998003 NA
4 season 1998 1998006 NA
5 season 1998 1998009 NA
6 season 1998 1998012 NA
Is there maybe some dependence used only for CPT that I should try to update?
Thank you!
day-time partititioning needs VPD as a predictor. Hence, if there are long gaps in VPD which are now filled by MDS algorithm, then also there are also long gaps in GPP_DT.
Hence, provide a method to estimate VPD in these long gaps from minimum temperature and dewpoint.
In very rare cases, the LRC fit is successful, but the inversion of the Hessian stops with an error.
Do catch this error and return an NA-result for this fit in LightResponseCurveFitter_fitLRC
.
Hello
We have been successfully running v 0.8-2 but recently we installed 1.0.0 on a new computer and are running into an error midway through our gapfill/partition script:
Start flux partitioning for variable NEE_NoUstar_f with temperature Tair_f.
Error in .self$sCalcPotRadiation(useSolartime.b = !isTRUE(debug.l$useLocaltime.b)) :
Need to set valid location information (sSetLocationInfo) before calling sCalcPotRadiation.
When running "sMRFluxPartition"
> EddyProc.C$sMRFluxPartition(Lat_deg.n=latitude, Long_deg.n=longitude, TimeZone_h.n=TimeZone, Suffix.s='NoUstar',
+ debug.l = list(useLocaltime.b = TRUE))
Start flux partitioning for variable NEE_NoUstar_f with temperature Tair_f.
Error in .self$sCalcPotRadiation(useSolartime.b = !isTRUE(debug.l$useLocaltime.b)) :
Need to set valid location information (sSetLocationInfo) before calling sCalcPotRadiation.
We noticed that we now need to run:
EddyProc.C$sSetLocationInfo(Lat_deg.n=latitude, Long_deg.n=longitude, TimeZone_h.n=TimeZone)
Which we have done but we still get this error. This code used to work so it seems we may be missing a flag/switch or some other setting to get this to run. We have tried to add/remove the debug flag which hasn't fixed it for us
Any help would be much appreciated. Thanks!
Currently, the scenario for which to keep all the intermediate results (argument uStarScenKeep
in sEddyProc_sGLFluxPartitionUStarScens
) defaults to character(0)
. This corresponds to the former convention of assigning no suffix to the default uStar scenario.
However, this changed with REddyProc Version to be reserved for the case without uStar filtering and using suffix "uStar"
for the case of uStarFiltering with unbootstrapped dataset.
adapt sEddyProc_sGLFluxPartitionUStarScens
default to use the first suffix, i.e. usually the "uStar"
scenario. User can explicitly supply "U50" to rather get the intermediates for the median boostraped uStar scenario.
The output of sEddyProc_sEstUstarThresholdDistribution
is a data.frame while the output of sEstUstarThreshold
is a list.
It makes it confusing now.
Actually, I would expect that if any of these two could be more complex and thus list would be needed though it is not. Maybe you could include in both methods argument fullOutput = FALSE producing just the expected data frame and list when set to TRUE. This way experienced users could explore more stats - if produced for $sEstUstarThresholdDistribution() - that are currently not available.
There are continuous support request on how to aggregate the uncertainty to daily and annual aggregates. Would be good to provide an example with the package.
I am trying to import from the AmeriFlux datasets (US-TW3 specifically) and I'm struggling with some of the columns. For example, the US-TW3 dataset doesn't include NEE. Is it possible to add an example similar to the AmeriFlux (US-TW3) data?
Thank you!
Currently sMDSGapFillAfterUstar
calls .self$sMDSGapFill
. If a user wants to modify the gapfilling, e.g. by including more environmental variables, then he has to recode the entire uStar-Filtering/Distribution procedure.
Implement an argument fGapFill
so that the user can provide its own modified gapfilling function to sMDSGapFillAfterUstar
.
This feature is used by clients, who do not need REddyProc. So it deserves its own package.
Currently the output-format written by fWriteDataframeToFile
is handled by fLoadTXTIntoDataframe
but hard to parse by other user code, maybe in other programming languages. This is largely due to keeping legacy compatibility.
Please, provide a format that is easier to parse.
Specifically
DateTime
or Date_Time
instead of Date Time
.I'm trying to process Ameriflux half-hourly flux data, for example for site US-CS1 and the file:
AMF_US-CS1_BASE_HH_2-5.csv
Which can be downloaded from Ameriflux. These data are indeed half-hourly and, though it's inconvenient, I used expand.grid()
and a left_join()
to guarantee that the data are on half-hourly steps (i.e., data are joined to a half-hourly time series by their date and hour of day).
require(REddyProc)
filename <- '~/Downloads/Ameriflux/extracts/AMF_US-CS1_BASE_HH_2-5.csv'
# Extract time of day, normalize field names
df0 <- read.csv(filename, skip = 2) %>%
select(date = TIMESTAMP_START, NEE = NEE_PI_F, Rg = SW_IN_1_1_1,
Tair = TA_1_1_1, VPD = VPD_PI_1_1_1, Ustar = USTAR_1_1_1) %>%
mutate(time = substr(date, 9, 12)) %>%
mutate(hour = as.integer(substr(time, 1, 2)) + as.integer(substr(time, 3, 4)) / 60) %>%
mutate(date = as.Date(substr(as.character(date), 1, 9), format = '%Y%m%d'))
# Create a continuous, half-hourly time series
df <- expand.grid(date = seq.Date(from = min(df0$date), to = max(df0$date), by = 'day'),
hour = seq(0, 23.5, 0.5)) %>%
# Join the original data to this series
left_join(df0, by = c('date', 'hour')) %>%
arrange(date, hour) %>%
# Extract meaningful date/time components for use in fConvertTimeToPosix()
mutate(year = format(date, '%Y'), month = format(date, '%m'), day = format(date, '%d'),
minute = ifelse(hour %% 1 == 0, 30, 0)) %>%
mutate(hour = floor(hour)) %>%
arrange(date, hour) %>%
mutate_at(.vars = vars(year, month, day,), as.integer) %>%
select(date, year, month, day, hour, minute, everything()) %>%
select(-time) %>%
# Fill in missing data with NA
mutate_all(.funs = function (x) ifelse(x == -9999, NA, x)) %>%
# Mask negative short-wave radiation values
mutate(Rg = ifelse(Rg < 0, NA, Rg)) %>%
fConvertTimeToPosix(
TFormat = 'YMDHM', Year = 'year', Month = 'month', Day = 'day', Hour = 'hour', Min = 'minute')
We can confirm the first few rows (at least) are half-hourly:
head(df)
> head(df)
DateTime date year month day hour minute NEE Rg Tair VPD Ustar
1 2018-01-01 00:30:00 17532 2018 1 1 0 30 NA NA NA NA NA
2 2018-01-01 00:00:00 17532 2018 1 1 0 0 NA NA NA NA NA
3 2018-01-01 01:30:00 17532 2018 1 1 1 30 NA NA NA NA NA
4 2018-01-01 01:00:00 17532 2018 1 1 1 0 NA NA NA NA NA
5 2018-01-01 02:30:00 17532 2018 1 1 2 30 NA NA NA NA NA
6 2018-01-01 02:00:00 17532 2018 1 1 2 0 NA NA NA NA NA
However, when I try to create a instance of sEddyProc
...
sEddyProc$new('US-CS1', df)
I get the error message:
Error in fCheckHHTimeSeries(Data[[ColPOSIXTime]], DTS = DTS, "sEddyProc.initialize") :
sEddyProc.initialize:::fCheckHHTimeSeries::: Time stamp is not equidistant (half-)hours in rows: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 2
Currently, the image plot of the fingerprint plot raises an error when there are infinite values in the data.
The infinite values should be replace by NA before plotting.
When implementing this feature, add unit tests and make the fingerprint work also outside of the reference class.
Currently fConvertTimeToPosix
relies on strptime
to return NA for implausible doy (day-of-year).
This causes a warning or worse on some architectures.
We should explicitly check for leapyears and compare each doy exceeding 365 or 366 respectively
before calling strptime
.
currently the user has to retrieve details, get the annual aggregation and set it again to the class.
Provide a method that does all this.
Currently daytime partitioning has to estimate temperature sensitivity, E0, and a prior for basal respiration, Rb, from night-time data. While the user can already specify E0, partitioning fails due to failing estimating prior of Rb from night-time data.
Supporting daytime partitioning in the absence of night-time records is important to process data from arctic stations where there is no night-time data.
Estimating the uncertainty of GPP with including the uncertainty in temperature sensitivity E0 requires a bootstrap of the data in each period, which consumes many resources. Currently, this is done for each uStar scenario in sEddyProc_sGLFluxPartitionUStarScens
. However, usually only the uncertainty is inspected for only one scenario.
The uncertainty estimate is used to determine, whether parameter are within bounds, but can be computed in a cheaper manner from the Hessian of the fit. Although this is a conservative estimate, that neglects uncertainty in E0, this is sufficient for the check of parameter bounds.
Omit the expensive uncertainty estimate for other uStar scenarios except one user-specified scenario.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.