This is the script for "Restricting the nonlinearity parameter in a commonly used soil greenhouse gas flux calculation scheme for more reliable flux estimates" published on PLOSone
| ############################################### | |
| ###### PAPER # 3 - Gas flux analysis ######### | |
| ############################################### | |
| # library(devtools) | |
| # install_bitbucket("ecoRoland/gasfluxes") | |
| library(gasfluxes) | |
| library(data.table) | |
| library(ggplot2) | |
| library(plot3D) | |
| library(plotly) | |
| # simulateflux <- function(t, V, A, c0 = 325, fmax = 5, nflux = 1, kappa = c(0.0001/36,exp(seq(-11,-3.5,0.32))), # fmax = 100/36, kappa was multiplied by 3600 | |
| simulateflux <- function(t, V, A, c0 = 325, fmax = 5, nflux = 25, kappa = c(0.0001/36,exp(seq(-11,-3.5,0.32))), # fmax = 100/36, kappa was multiplied by 3600 | |
| nMC , GC.sd, g.factor = 4, f.detect, | |
| seed) { | |
| if (!missing(seed)) set.seed(seed) | |
| ## model input data ## | |
| h <- V/A | |
| # f0 <- signif(seq(0,fmax,length.out = nflux),3) # all fluxes to be simulated from 0 to fmax; | |
| # f0 <- signif(exp(seq(0,log(fmax+1),length.out = nflux)),5)-1 # all fluxes to be simulated from 0 to fmax; | |
| f0 <- signif(c(0,exp(seq(-5,log(fmax),length.out=nflux-1))),3) # exponential f0 scale to have more data on small fluxes | |
| # if (missing(f.detect)) { | |
| # if (length(t) == 4L) { | |
| # | |
| # # f.detect <- 13.20 * t[4]^(-0.9973) * c0 * GC.sd/c0 * A | |
| # # f.detect <- 0.036 # calculated with GC.sd = 5 | |
| # f.detect <- 0.0808 # 0.0808 nmol N2O / m2 / s for GD.sd = 3 (always HMR) | |
| # f.detect <- 0.08000 # nmol N2O / m2 / s for GD.sd = 3 equilibrium | |
| # | |
| # } else { | |
| #simulate f.detect | |
| f.detect <- 0 # starting value for f.detect | |
| C00 <- c0 * 273.15 / 22.4 / (273.15 + 15) * 1000 #nmol N2O / m^3 | |
| sdGC <- GC.sd * 273.15 / 22.4 / (273.15 + 15) * 1000 #from ppb GC.SD to nmol N / m^3 | |
| # t = c(seq(0,01,length.out=4))*60 # c(0,1900,3000,3600) | |
| for (i in 1:2) { | |
| set.seed(24) | |
| # sim <- data.frame(t = seq(0, 1, length.out = length(t)), C = rnorm(4e3, mean = C00, sd = sdGC), | |
| # id = rep(1:1e3, each = length(t)), A = A, V = V) | |
| sim <- data.frame(t = t, C = rnorm(4e3, mean = C00, sd = sdGC), | |
| id = rep(1:1e3, each = length(t)), A = A, V = V) | |
| simflux <- gasfluxes(sim, .id = "id", .times = "t", methods = c("HMR", "linear"), k_HMR = log(1.5/3600), plot = FALSE) | |
| simflux[, f0 := linear.f0] # use linear estimate as default for this method | |
| # simflux[, t.inf := (HMR.phi-linear.C0)/linear.f0/h ] # sec mean(t) | |
| # simflux[, kappa.max := abs(linear.f0/f.detect/max(t))] | |
| # simflux[, kappa.max := linear.f0/f.detect / (t.inf - sum(t)) ] # tfactor (t.inf - sum(t))/3600 | |
| simflux[ !is.na(HMR.f0), f0 := HMR.f0] # | |
| n.HMR <- simflux[ !is.na(HMR.f0), .N] # | |
| f.detect <- simflux[, quantile(f0, 0.975)] # 0.0808 nmol N2O / m2 / s for GD.sd = 3 (always HMR) | |
| write.table(t(c(as.numeric(f.detect),GC.sd,n.HMR)),file="fdetect-out.csv",append=T,col.names = F, sep=";") | |
| } | |
| # selectfluxes(res, "RH2017", f.detect = f.detect) | |
| # res[method == "HMR", .N] #9 | |
| # | |
| # stop("Detection limit can only be estimated for flux measurements with three or four time points. Please specify the f.detect parameter.") | |
| # } | |
| # } | |
| # Ct.test <- (c0 - fmax/(-tail(kappa, 1) * h)) + fmax * exp(-tail(kappa, 1)*t)/(-tail(kappa, 1)*h) # in c0 in ppb | |
| Ct.test <- (C00 - fmax/(-tail(kappa, 1) * h)) + fmax * exp(-tail(kappa, 1)*t)/(-tail(kappa, 1)*h) # in c0 in ppb | |
| # neftel.rsq <- signif(summary(lm(Ct.test ~t))$r.squared,5) # calculation of R.lim (minimal R.squared for a given system) | |
| # prepare input for flux calculation for all simulated fluxes | |
| df <- data.table(V=V, A=A, time=t) | |
| df <- df[rep(seq_len(nrow(df)), nMC)] | |
| df[, mc := rep(seq_len(nMC), each = length(t))] | |
| df <- df[rep(seq_len(nrow(df)), length(f0) * length(kappa))] | |
| combs <- CJ(f0 = f0, kappa = kappa) | |
| df[, c("f0", "kappa") := combs[rep(seq_len(nrow(combs)), each = length(t) * nMC)]] | |
| df[, phi := c0 - f0/(-kappa*h)] | |
| df[, C := phi + f0 * exp(-kappa*t)/(-kappa*h) + rnorm(.N, sd = GC.sd)] | |
| #flux calculation | |
| message("Calculating ", nrow(df) / length(t), " fluxes. Please be patient... jaah ganz easy!") | |
| res <- gasfluxes(df, c("mc", "f0","kappa", "phi"),methods = c("linear", "robust linear", "HMR"), gfactor = g.factor, plot = FALSE, | |
| verbose = FALSE,k_HMR= log(1.5/3600)) | |
| # res[, rsq := df[, cor(C, time)^2, by = c("mc", "kappa", "phi")][["V1"]]] #have gasfluxes return rsq instead | |
| ############################### | |
| ### decision making section ### | |
| ############################### | |
| # RF2011new | |
| res <- selectfluxes(res, select = "RF2011new") | |
| setnames(res,"flux","RF2011")#[] | |
| # ## dynamic R.squared method, dyn.rsq ## | |
| # f.detect.rh <- res[f0 == 0, quantile(linear.f0, 0.95)] # detection limit is 95 % quantile of all linear zero fluxes, correct?????? | |
| # res[, rsq.lim := ((3*f.detect.rh)^2 + neftel.rsq*linear.f0^2) / ((3*f.detect.rh)^2 + linear.f0^2)] # formula that calculates the R² threshold allowing for nonlinear fluxes; reformulated dyn.rsq definition by Roland | |
| # res[, range(rsq.lim)] | |
| # res[, dyn.rsq := linear.f0] | |
| # res[rsq > rsq.lim & !is.na(HMR.f0), dyn.rsq := HMR.f0] # use HMR estimate if HMR was fitted and method requirements are met | |
| # | |
| ### implemented in the selectflux function ### | |
| # res <- selectfluxes(res, "RH2017", f.detect = f.detect, beta = 3, R.lim = neftel.rsq) | |
| # res[, dyn.Rsq := flux] | |
| # res <- selectfluxes(res, "RH2017", f.detect = f.detect, | |
| # dyn.Rsquared = function(f.lin, R.lim, f.detect) ((30 * f.detect)^2 + R.lim * f.lin^2)/((30 * f.detect)^2 + f.lin^2) | |
| # ) | |
| # ## minimum detectable flux scheme, f.det using f.detect ## | |
| # res[, f.det := linear.f0] # use linear estimate as default for this method | |
| # res[linear.f0 > 2*f.detect & !is.na(HMR.f0), f.det := HMR.f0] | |
| # | |
| ## g-factor scheme, g.fact ## | |
| res[, g.fact := robust.linear.f0] # use robust linear estimate as default for this method | |
| res[abs(HMR.f0/linear.f0)<g.factor & !is.na(HMR.f0), g.fact := HMR.f0] #note the abs!!! | |
| res[is.na(g.fact), g.fact := linear.f0] # use linear if robust and HMR dont work | |
| ## p-value scheme, p.value ## | |
| # res[, p.value := linear.f0] # use linear estimate as default for this method | |
| # res[ HMR.f0.p < linear.f0.p & !is.na(HMR.f0), p.value := HMR.f0] #note the abs!!! | |
| # | |
| ## fix Rsq ## | |
| # res[, Rsq.fix := linear.f0] # use linear estimate as default for this method | |
| # res[ 0.8 < (linear.r)^2 & !is.na(HMR.f0), Rsq.fix := HMR.f0] #note the abs!!! | |
| ## fix kappa ## | |
| # res[, kappa.fix := linear.f0] # use linear estimate as default for this method | |
| # res[ 5/3600 > HMR.kappa & !is.na(HMR.f0), kappa.fix := HMR.f0] #note the abs!!! | |
| # | |
| ## dynamic kappa scheme, kappa.max ## | |
| res[, kappa.max := robust.linear.f0] # use robust linear estimate as default for this method | |
| res[, kappa.Max := linear.f0/f.detect/max(t)] | |
| res[ HMR.kappa < kappa.Max & !is.na(HMR.f0), kappa.max := HMR.f0] # | |
| res[is.na(kappa.max), kappa.max := linear.f0] # use linear if robust and HMR dont work | |
| # ## dynamic kappa scheme, kappa.max + AIC ## | |
| # res[, kappa.AIC := linear.f0] # use linear estimate as default for this method | |
| # res[, kappa.AIC := linear.f0/f.detect/max(t)] | |
| # res[ HMR.kappa < kappa.AIC & !is.na(HMR.f0), kappa.AIC := HMR.f0] # | |
| ## dynamic kappa scheme PLUS, dyn.kpp ## | |
| # res[, dyn.kpp := linear.f0] # use linear estimate as default for this method res[, t.inf := (HMR.phi-linear.C0)/linear.f0/h ] # sec mean(t) | |
| # res[, t.inf := (HMR.phi-linear.C0)* h / linear.f0] # sec corrected version | |
| # # res[, t.inf := (HMR.phi-linear.C0)/linear.f0/h ] # sec by h instead of * h | |
| # res[, kpp.max := linear.f0 / f.detect / (t.inf- max(t)) ] # ## tfactor (t.inf - sum(t))/3600 | |
| # # res[, kpp.max := 1 / (t.inf) ] # tfactor (t.inf - sum(t))/3600 | |
| # # res[, kpp.max := abs(linear.f0/f.detect) / (var(t) / max(t)) ] # variance into t.factor introduced | |
| # res[ HMR.kappa < kpp.max & !is.na(HMR.f0), dyn.kpp := HMR.f0] # | |
| ## Rolands new scheme, fuss.2017 ## | |
| res[, RF2017 := robust.linear.f0] # use linear estimate as default for this method | |
| res[HMR.kappa < (20/3600) & !is.na(HMR.f0), RF2017 := HMR.f0] # Rolands new criteria explained by email, kappa smaller 5 | |
| res[is.na(RF2017), RF2017 := linear.f0] # use linear if robust and HMR dont work | |
| ## Thus I used HMR if it could be fit and kappa was below 20 h^-1 (in practice this resulted in kappa < 5 for all fluxes accepted as non-linear). ## | |
| ## AIC scheme ## | |
| res[, AIC := robust.linear.f0] # use linear estimate as default for this method | |
| res[ HMR.AIC < linear.AIC & !is.na(HMR.f0), AIC := HMR.f0] # & HMR.kappa < 20/3600 ; Rolands new criteria explained by email, kappa smaller 5 | |
| res[is.na(AIC), AIC := linear.f0] # use linear if robust and HMR dont work | |
| ## HMR-linear combination scheme, HMR.f0 ## | |
| res[is.na(HMR.f0), HMR.f0 := linear.f0] # Use linear for HMR if no HMR could be fitted | |
| res[is.na(robust.linear.f0), robust.linear.f0 := linear.f0] | |
| class(res) <- c("gasfluxes.simu", class(res)) | |
| res[] #return value | |
| } | |
| ############################### | |
| ### start actual simulation ### | |
| # simus <- simulateflux(t = c(seq(0,36,length.out=4))/60, V = 0.0085, A = 0.0707, nMC = 10, seed = 42) | |
| # system.time(simus <- simulateflux(t = c(seq(0,60,length.out=4))*60, V = 0.015, A = 0.07, GC.sd = 3, nMC = 25, seed = 88)) # fibl system, takes 3.5 min for nMC = 40 | |
| system.time(simus <- simulateflux(t = c(seq(0,36,length.out=4))*60, V = 0.014, A = 0.07, GC.sd = 3, nMC = 100, seed = 12)) # Nora system, takes 3.5 min for nMC = 40 | |
| names(simus)[30] <- "dyn.Rsq" | |
| # simus.55.backup <-simus.45.backup | |
| # simus.40.backup <- simus | |
| # simu.rf2011 <- selectfluxes(simus, "RF2011new") | |
| # simus$g.fact <- simu.rf2011$flux | |
| # dyn.Rsquared = function(f.lin, R.lim, f.detect) ((30 * f.detect)^2 + R.lim * f.lin^2)/((30 * f.detect)^2 + f.lin^2) | |
| # simu.dynRsqpls <- selectfluxes(simus, "RH2017", f.detect = f.detect) | |
| #relatively slow: | |
| # system.time(simus <- simulateflux(t = c(seq(0,36,12))/60, V = 0.0085, A = 0.0707, nMC = 100, seed = 42)) | |
| #on my (fast) system: 400 seconds | |
| # simus[HMR.f0 > 3.3, HMR.f0 := 3.3] # ??? what was that good for? ... limit HMR.f0 estimate on 3.3?? | |
| summary.gasfluxes.simu <- function(x) { | |
| x <- melt(x, id.vars = c("mc", "f0", "kappa", "phi"), | |
| # measure.vars = c("linear.f0","robust.linear.f0","HMR.f0","g.fact","kappa.max","AIC"), | |
| measure.vars = c("linear.f0","robust.linear.f0","HMR.f0","AIC","g.fact","kappa.max","RF2011","RF2017"), # | |
| # measure.vars = c("linear.f0","kappa.max","HMR.f0"), | |
| variable.name = "method", value.name = "flux.estimate") | |
| x <- x[, .(median = median(flux.estimate), | |
| mean = mean(flux.estimate), | |
| sd = sd(flux.estimate,na.rm = T), | |
| lwr = quantile(flux.estimate, 0.05,na.rm=T), | |
| upr = quantile(flux.estimate, 0.95,na.rm=T), | |
| var = var(flux.estimate,na.rm = T)), by = c("f0", "kappa", "phi", "method")] | |
| x[, IQR := upr - lwr] | |
| x[, bias:= mean - f0] | |
| x[, MSE := bias^2 + var] | |
| class(x) <- c("gasfluxes.simu.summary", class(x)) | |
| x[] | |
| } | |
| contour.fix <- function(df) { | |
| colname <- names(df) | |
| names(df) <- c("x", "y", "z") | |
| Range <- as.data.frame(sapply(df, range)) | |
| Range[1,1] <- 1E-6 | |
| Dim <- as.data.frame(t(sapply(df, function(x) length(unique(x))))) | |
| arb_z = Range$z[1] - diff(Range$z)/Dim$x | |
| df2 <- rbind(df, | |
| expand.grid(x = c(exp(log(Range$x[1]) - diff(log(Range$x))/Dim$x), exp(log(Range$x[2]) + diff(log(Range$x))/Dim$x)), | |
| y = unique(df$y), z = arb_z), # seq(Range$y[1], Range$y[2], length = Dim$y) | |
| expand.grid(x = unique(df$x), # seq(Range$x[1], Range$x[2], length = Dim$x) | |
| y = c(exp(log(Range$y[1]) - diff(log(Range$y))/Dim$y), exp(log(Range$y[2]) + diff(log(Range$y))/Dim$y)), z = arb_z)) | |
| # g <- ggplot(df2, aes(x, y, z = z)) + labs(x = colname[1], y = colname[2], fill = colname[3]) + | |
| # stat_contour(geom="polygon", aes(fill=..level..),bins =50) + | |
| # coord_cartesian(xlim=c(Range$x), ylim=c(Range$y), expand = F) + | |
| # scale_y_log10(breaks = c(1E-2,1E-3,1E-4,1E-5)) + | |
| # theme_bw() | |
| names(df2) <- c("f0", "kappa", "stat.x") | |
| return(df2) | |
| } | |
| simus.sum <- summary(simus) | |
| # simus.sumus <- summary.gasfluxes.simu(simus) | |
| # head(simus.sum) | |
| # head(simus.sumus) | |
| # median(simus$dyn.Rsq) | |
| # simus.sum[IQR > 0.5, IQR := 0.5] # Use linear for HMR if no HMR could be fitted | |
| # simus.sum[, IQR := log(IQR)] # use logarithmic IQR values | |
| # simus.sum[, var := log(var)] # use logarithmic variance values | |
| # simus.sum[, MSE := log(MSE)] # use logarithmic MSE values | |
| plot.gasfluxes.simu.mean <- function(x, type = "contour", stat = "median") { | |
| if (type == "contour") { | |
| simus.plot <- simus.sum[,contour.fix(.(f0,kappa,median)),method] | |
| p <- ggplot(simus.plot, aes_(x = ~f0, y = ~kappa, z = ~stat.x)) + | |
| stat_contour(aes(fill = ..level..), geom = "polygon", bins = 128) + | |
| facet_wrap(~ method) + | |
| scale_fill_distiller(name = sprintf("flux estimate\n(%s)", stat), palette = "BrBG") + # BRBG "YlGnBu" RdYlGn" | |
| scale_y_log10(breaks = c(1E-2,1E-3,1E-4,1E-5)) + | |
| coord_cartesian(xlim = range(simus.sum$f0), ylim= range(simus.sum$kappa), expand = F) + | |
| theme_bw() | |
| return(print(p)) | |
| } | |
| } | |
| plot.gasfluxes.simu.mean(simus.sum ) | |
| plot.gasfluxes.simu.iqr <- function(x, type = "IQR.plot", stat = "IQ90") { | |
| if (type == "IQR.plot") { | |
| simus.plot <- simus.sum[,contour.fix(.(f0,kappa,log(IQR))),method] | |
| p <- ggplot(simus.plot, aes_(x = ~f0, y = ~kappa, z = ~stat.x)) + | |
| stat_contour(aes(fill = ..level..), geom = "polygon", bins = 128) + | |
| facet_wrap(~ method) + | |
| scale_fill_distiller(name = sprintf("estimated\n log(%s)", stat), palette = "RdYlGn") + # BrBG "YlGnBu" RdYlGn" | |
| scale_y_log10(breaks = c(1E-2,1E-3,1E-4,1E-5)) + | |
| coord_cartesian(xlim = range(simus.sum$f0), ylim= range(simus.sum$kappa), expand = F) + | |
| theme_bw() | |
| # annotate("text", x = 2, y=1E-5, label = c("blablab")) | |
| return(print(p)) | |
| } | |
| } #more plot types can be implemented | |
| plot.gasfluxes.simu.iqr(simus.sum) | |
| plot.gasfluxes.simu.var <- function(x, type = "var.plot", stat = "var") { | |
| if (type == "var.plot") { | |
| simus.plot <- simus.sum[,contour.fix(.(f0,kappa,log(var))),method] | |
| p <- ggplot(simus.plot, aes_(x = ~f0, y = ~kappa, z = ~stat.x)) + | |
| stat_contour(aes(fill = ..level..), geom = "polygon", bins = 128) + | |
| facet_wrap(~ method) + | |
| scale_fill_distiller(name = sprintf("estimated\n log(%s)", stat), palette = "RdYlGn") + # BrBG "YlGnBu" RdYlGn" | |
| scale_y_log10(breaks = c(1E-2,1E-3,1E-4,1E-5)) + | |
| coord_cartesian(xlim = range(simus.sum$f0), ylim= range(simus.sum$kappa), expand = F) + | |
| theme_bw() | |
| # annotate("text", x = 2, y=1E-5, label = c("blablab")) | |
| return(print(p)) | |
| } | |
| } #more plot types can be implemented | |
| plot.gasfluxes.simu.var(simus.sum) | |
| plot.gasfluxes.simu.var <- function(x, type = "var.plot", stat = "MSE") { | |
| if (type == "var.plot") { | |
| simus.plot <- simus.sum[,contour.fix(.(f0,kappa,log(MSE))),method] | |
| p <- ggplot(simus.plot, aes_(x = ~f0, y = ~kappa, z = ~stat.x)) + | |
| stat_contour(aes(fill = ..level..), geom = "polygon", bins = 128) + | |
| facet_wrap(~ method) + | |
| scale_fill_distiller(name = sprintf("estimated\n log(%s)", stat), palette = "RdYlGn") + # BrBG "YlGnBu" RdYlGn" | |
| scale_y_log10(breaks = c(1E-2,1E-3,1E-4,1E-5)) + | |
| coord_cartesian(xlim = range(simus.sum$f0), ylim= range(simus.sum$kappa), expand = F) + | |
| theme_bw() | |
| # annotate("text", x = 2, y=1E-5, label = c("blablab")) | |
| return(print(p)) | |
| } | |
| } #more plot types can be implemented | |
| plot.gasfluxes.simu.var(simus.sum) | |
| simus.sum[,mean(exp(IQR)),method] # average IQR for each method | |
| simus.sum[f0 == 0, .(quantile(mean,0.975)), method] # detection limit for each method | |
| ### ---- | |
| # plot(nora.dataset$flux,log(nora.dataset$kappa), xlim = c(0,2.78), ylim = c(1E-5,1E-2),col="black") | |
| # plot.gasfluxes.hist(simus.sum) | |
| # simus.sum[simus.sum$method=="dyn.Rsq" & simus.sum$f0 == 0.333 ] | |
| simus.sum.backup <- simus.sum | |
| # | |
| # simus.sum.tmax13 <- simus.sum | |
| # | |
| # | |
| # simus.60 <- simus.sum.tmax60[method=="kappa.max"] | |
| # simus.60[,method:="t.max = 60 min"] | |
| # simus.45 <- simus.sum.tmax45[method=="kappa.max"] | |
| # simus.45[,method:="t.max = 45 min"] | |
| # simus.30 <- simus.sum.tmax30[method=="kappa.max"] | |
| # simus.30[,method:="t.max = 30 min"] | |
| # simus.12 <- simus.sum.tmax12[method=="kappa.max"] | |
| # simus.12[,method:="t.max = 12 min"] | |
| # simus.75 <- simus.sum.tmax75[method=="kappa.max"] | |
| # simus.75[,method:="t.max = 75 min"] | |
| # simus.90 <- simus.sum.tmax90[method=="kappa.max"] | |
| # simus.90[,method:="t.max = 90 min"] | |
| # simus.GCsd <- rbind(simus.15,simus.09,simus.06,simus.03,simus.01) | |
| # simus.GCsd <- rbind(simus.13,simus.30,simus.45,simus.60,simus.75,simus.90) | |
| # simus.sum <- simus.GCsd | |
| # plot.gasfluxes.simu.iqr(simus.sum) | |
| ########################################### | |
| #### density frequence comparision ######## | |
| ########################################### | |
| # load a GHG flux dataset... | |
| # lysi.gfg <- read.table("gfg_Lysi.csv",sep=",",header=T) | |
| ## add fibl data | |
| # fibl.raw <- read.table("/home/roman/PhD/data/GasFluxGames/Fibl/N2O_Rohdaten.csv",header=T, sep = ",") | |
| # fibl.raw <- read.table("/home/hueppir/nashome/Gasflux-Simulations/GasFluxGames/Fibl/N2O_Rohdaten.csv",header=T, sep = ",") | |
| fibl.raw <- fread("/home/hueppir/naslocal/Gasflux-Simulations/GasFluxGames/Fibl/N2O_Rohdaten.csv",header=T, sep = ",") | |
| fibl.raw$time <- fibl.raw$time * 3600 # we use seconds instead of hours | |
| # art.raw <- fread("/home/hueppir/PhD/data/Reckenholz2014/fluxes-n2o/n2o-input-complete.csv",header=T,sep = ",") | |
| # art.raw$V <- art.raw$V / 1000 | |
| # fibl.raw <- art.raw | |
| fibl.raw$C <- fibl.raw$C * 1E3 / 22.4136 * (273.15/(15 + 273.15)) # n mol N2O/m^3 ;22.4L/mol;273K, | |
| # flux.ghg <- gasfluxes(fibl.raw, methods = c("linear", "robust linear", "HMR"), k_HMR = log(1.5/3600),plot = T) # fluxes in nmol-N2O/m²/s | |
| flux.ghg <- gasfluxes(na.omit(fibl.raw), | |
| methods = c("linear", "robust linear", "HMR"), | |
| .id = c("ID", "V", "A"), | |
| .V = "V", .A = "A", .times = "time", .C = "C", k_HMR = log(1.5/3600),plot = F) # fluxes in nmol-N2O/m²/s | |
| # plot_ly(flux.ghg,x = flux.ghg$linear.C0, y = flux.ghg$HMR.phi) | |
| # plot_ly(data=flux.ghg, x = flux.ghg$HMR.kappa, y = flux.ghg$t.phi ) | |
| # write.table(flux.ghg, "/home/roman/PhD/data/GasFluxGames/Fibl/fibl-n2o-fluxes.csv", sep = ";", row.names = F) | |
| write.table(flux.ghg, "/home/hueppir/naslocal/Gasflux-Simulations/GasFluxGames/Fibl/fibl-n2o-fluxes_detailed.csv", sep = ";", row.names = F) | |
| # write.table(flux.ghg, "/home/hueppir/naslocal/Gasflux-Simulations/GasFluxGames/Reckenholz/art-n2o-fluxes_detailed.csv", sep = ";", row.names = F) | |
| flux.ghg.backup <- flux.ghg | |
| flux.ghg <- fread("/home/hueppir/naslocal/Gasflux-Simulations/GasFluxGames/Fibl/fibl-n2o-fluxes_detailed.csv") | |
| # always HMR if fitted | |
| flux.ghg[is.na(HMR.f0), HMR.f0 := linear.f0] # Use linear for HMR if no HMR could be fitted | |
| flux.ghg[is.na(HMR.f0), HMR.f0 := robust.linear.f0] | |
| flux.ghg[is.na(robust.linear.f0) & is.na(HMR.f0.se),HMR.f0 := linear.f0] | |
| flux.ghg[HMR.f0 > 5000, HMR.f0 := linear.f0] ## avoid the very large HMR fluxes | |
| ### kappa.max method | |
| f.detect <- 0.069 # from calculation 0.028 | |
| t <- c(seq(0,36,length.out=4))*60 | |
| flux.ghg[, kappa.max := robust.linear.f0] # use linear estimate as default for this method | |
| flux.ghg[is.na(kappa.max), kappa.max := linear.f0] # use linear estimate as default for this method | |
| flux.ghg[is.na(robust.linear.f0) & is.na(HMR.f0.se),kappa.max := linear.f0] | |
| flux.ghg[, kappa.dyn := abs(linear.f0/f.detect/max(t))] | |
| flux.ghg[ HMR.kappa < kappa.dyn & !is.na(HMR.f0.se), kappa.max := HMR.f0] # | |
| # flux.ghg[HMR.kappa < kappa.max & !is.na(HMR.f0.se), .N] | |
| # flux.ghg[,mean(dyn.kappa,na.rm = T)] | |
| ### AIC method | |
| flux.ghg[, AIC := robust.linear.f0] # use linear estimate as default for this method | |
| flux.ghg[ HMR.AIC < linear.AIC & !is.na(HMR.f0.se), AIC := HMR.f0] # | |
| flux.ghg[is.na(robust.linear.f0.se) & is.na(HMR.f0.se),AIC := linear.f0] | |
| # flux.ghg[ HMR.AIC < linear.AIC & !is.na(HMR.f0.se), .N] | |
| # flux.ghg[,mean(AIC,na.rm = T)] | |
| ### g.fact | |
| flux.ghg[, g.fact := robust.linear.f0] # use linear estimate as default for this method | |
| flux.ghg[HMR.f0/robust.linear.f0 < 4 & !is.na(HMR.f0.se), g.fact := HMR.f0] # Rolands new criteria explained by email, kappa smaller 5 | |
| flux.ghg[is.na(g.fact), g.fact := linear.f0] # Use linear for HMR if no HMR could be fitted | |
| # flux.ghg[HMR.f0/robust.linear.f0 < 4 & !is.na(HMR.f0.se) , .N] | |
| # flux.ghg[,mean(g.fact,na.rm = T)] | |
| ## Rolands new scheme, fuss.2017 ## | |
| max.kp <- 20 | |
| flux.ghg[, RF2017 := robust.linear.f0] # use linear estimate as default for this method | |
| flux.ghg[HMR.kappa < (max.kp/3600) & !is.na(HMR.f0.se), RF2017 := HMR.f0] # Rolands new criteria explained by email, kappa smaller 5 | |
| flux.ghg[is.na(RF2017), RF2017 := linear.f0] # Use linear for HMR if no HMR could be fitted | |
| # flux.ghg[HMR.kappa < max.kp/3600 & !is.na(HMR.f0.se), .N] | |
| # flux.ghg[,mean(fuss.2017,na.rm = T)] | |
| ### RF2011new; bad p-value criteria | |
| selectfluxes(flux.ghg, select = "RF2011new") | |
| setnames(flux.ghg,"flux","RF2011")[] | |
| ### Robust linear use linear if only 3 data points | |
| flux.ghg[is.na(robust.linear.f0), robust.linear.f0 := linear.f0] # Use linear for HMR if no HMR could be fitted | |
| # flux.ghg$robust.linear.f0[is.na(flux.ghg$robust.linear.f0)] <- flux.ghg$linear.f0[is.na(flux.ghg$robust.linear.f0)] | |
| # flux.ghg[method == "HMR", .N] | |
| # flux.ghg[,mean(flux,na.rm = T)] | |
| # method.tab$mean[3] <- flux.ghg[,mean(HMR.f0,na.rm = T)] | |
| # flux.ghg[!is.na(HMR.AIC), .N] | |
| # fibl.melt <- melt(flux.ghg, id.vars = c("ID"), | |
| # measure.vars = c("linear.f0", "HMR.f0","robust.linear.f0"), | |
| # variable.name = "method", value.name = "flux.estimate") | |
| # fibl.results <- fibl.melt[,.(mean.flux=mean(flux.estimate,na.rm = T)),method] | |
| # flux.ghg <- flux.ghg.backup | |
| # flux.ghg <- fread("/home/hueppir/nashome/Gasflux-Simulations/GasFluxGames/Reckenholz/art-n2o-fluxes_detailed.csv") | |
| # selectfluxes(flux.ghg, select = "RH2017", R.lim = 0.6, f.detect = 0.067) | |
| # flux.ghg <- selectfluxes(flux.ghg, select = "RH2017", R.lim = 0.352, f.detect = 0.39) | |
| # method.tab <- data.table(methods = c("linear","robust linear","HMR","AIC","g-factor","kappa.max","RF2011","RF2017")) | |
| # matrix(NA, nrow = 8, ncol = 6) | |
| # row.names(method.tab) <- c("linear","robust linear","HMR","AIC","g-factor","kappa.max","RF2011","RF2017") | |
| # names(method.tab) <- c("mean","n HMR","n linear","nonlinear/linear","fdet","MSE") | |
| # method.tab$mean <- NA | |
| # method.tab$mean[1] <- flux.ghg[,mean(linear.f0)] | |
| # method.tab$mean[2] <- flux.ghg[,mean(robust.linear.f0,na.rm=T)] | |
| # flux.ghg[, .(mean = mean)] | |
| # flux.ghg[method == "HMR", .N] | |
| # flux.ghg[,mean(flux,na.rm = T)] | |
| # HMR method, robust.linear alternatively | |
| # method.tab$nonlinear.linear <- NA | |
| # method.tab$bias <- NA | |
| # method.tab$IQ90 <- NA | |
| # method.tab$MSE <- NA | |
| # flux.ghg[is.na(robust.linear.f0),.N] | |
| # flux.ghg[is.na(robust.linear.f0),robust.linear.f0 := linear.f0] | |
| # method.tab$bias <- flux.ghg[,mean(robust.linear.f0,na.rm = T)] | |
| # h <- 0.2 | |
| # f.detect <- 0.04 | |
| # ### kpp.max (dyn.kappa) method that is more restrictive and uses t.meas equal to t.phi-t.meas | |
| # flux.ghg[, dyn.kappa.plus := robust.linear.f0] # use linear estimate as default for this method | |
| # flux.ghg[is.na(dyn.kappa.plus), dyn.kappa.plus := linear.f0] # use linear estimate as default for this method | |
| # flux.ghg[, t.inf := (HMR.phi-linear.C0)/linear.f0/h ] # sec mean(t) | |
| # flux.ghg[, kpp.max := linear.f0/f.detect / (t.inf - sum(t)) ] # tfactor (t.inf - sum(t))/3600 | |
| # # res[, kpp.max := abs(linear.f0/f.detect) / (var(t) / max(t)) ] # variance into t.factor introduced | |
| # flux.ghg[ HMR.kappa < kpp.max & !is.na(HMR.f0), dyn.kappa.plus := HMR.f0] # | |
| # | |
| # flux.ghg[HMR.kappa < kpp.max, .N] | |
| # flux.ghg[,mean(dyn.kappa.plus,na.rm = T)] | |
| # plot_ly(data=flux.ghg, x = flux.ghg$HMR.kappa, y = flux.ghg$HMR.phi ) | |
| geo_mean <- function(data) { | |
| min_data <- min(data,na.rm = T) | |
| log_data <- log(data-min_data) | |
| gm <- exp(mean(log_data[is.finite(log_data)]))+min_data | |
| return(gm) | |
| } | |
| # ### RF2011old; wrong p-values | |
| # flux.ghg.orig <- gasfluxes(fibl.raw, methods = c("linear", "robust linear", "original HMR"), k_HMR = log(1.5/3600),plot = F) # fluxes in nmol-N2O/m²/s | |
| # write.table(flux.ghg.orig, "/home/roman/PhD/data/GasFluxGames/Fibl/fibl-n2o-originalHMR-fluxes.csv", sep = ";", row.names = F) | |
| # selectfluxes(flux.ghg.orig, select = "RF2011") | |
| # flux.ghg.orig[method == "original HMR", .N] | |
| # flux.ghg.orig[,mean(flux,na.rm = T)] | |
| flux.ghg.DT <- flux.ghg | |
| write.table(flux.ghg, "/home/hueppir/naslocal/Gasflux-Simulations/GasFluxGames/Fibl/fibl-n2o-results_detailed2.csv", sep = ";", row.names = F) | |
| flux.ghg <- read.table("/home/hueppir/naslocal/Gasflux-Simulations/GasFluxGames/Fibl/fibl-n2o-results_detailed2.csv",sep=";",header = T) | |
| # flux.ghg <- read.table("L:/Gasflux-Simulations/GasFluxGames/Fibl/fibl-n2o-results_detailed2.csv",sep=";",header = T) | |
| # flux.ghg[,..cols] | |
| hist.eval <- function(data, fcs.method, #fcs.parameter = "IQR", | |
| fmax = 5, nflux = 25, kappa = c(0.0001/36,exp(seq(-11,-3.5,0.32))), | |
| f0 = signif(c(0,exp(seq(-5,log(fmax),length.out=nflux-1))),3)) { | |
| # substitute(x) | |
| flux.ghg <- as.data.frame(data) | |
| # fcs.method <- substitute(kappa.max) | |
| # flux.ghg <- read.table("F:/SAE/Matti/Nora/Nora-fluxes-n2o_all.csv",sep = ";",header = T) | |
| # # flux.ghg <- read.table("/home/roman/SAE/Matti/Nora/Nora-fluxes-n2o_all.csv",sep = ";",header = T) | |
| # fmax <- 5 # 100/36 | |
| # nflux <- 25 | |
| # kappa <- c(0.0001/36,exp(seq(-11,-3.5,0.32))) | |
| # f0 <- signif(c(0,exp(seq(-5,log(fmax),length.out=nflux-1))),3) # exponential f0 scale to have more data on small fluxes | |
| # kappa.na <- 0.0001 | |
| # kappa.na <- 0.0001/36 # = min(kappa) 1E-6 | |
| # flux.ghg$flux.hist <- flux.ghg$rsq.rh | |
| # nora.dataset <- data.frame(flux = flux.ghg$linear.f0, kappa= flux.ghg$HMR.kappa) ## prepare dataset for 3d plot combination | |
| # nora.dataset$kappa[is.na(nora.dataset$kappa)] <- 0.01 # NA flux values are 0 | |
| # simus[ ,nora.kappa:=flux.ghg$HMR.kappa , method] | |
| # mean.flux.dataset <- mean(flux.ghg$flux,na.rm=T) | |
| ### reframe measured data into model space ### | |
| flux.ghg$flux.hist <- flux.ghg[,fcs.method] # flux.ghg$linear.f0 ## insert actual method results here!! | |
| flux.ghg$flux.hist[is.na(flux.ghg$flux.hist)] <- 0 # NA flux values are 0 | |
| flux.ghg$flux.hist[flux.ghg$flux.hist < 0] <- abs(flux.ghg$flux.hist[flux.ghg$flux.hist < 0]) # transform negative to positive fluxes | |
| flux.ghg$flux.hist[flux.ghg$flux.hist>fmax] <- fmax ## 100 # fmax 2.8 % of Maikes fluxes are higher than that | |
| flux.ghg$kappa.hist <- flux.ghg$HMR.kappa | |
| kappa.hist <- kappa | |
| flux.ghg$kappa.hist[flux.ghg$HMR.diagnostics=="singulärer Gradient"] <- min(kappa.hist) # lower end of simulation "singular gradient" | |
| flux.ghg$kappa.hist[is.na(flux.ghg$kappa.hist)] <- max(kappa.hist)*0.95 # higher end of simulation | |
| # flux.ghg$kappa.hist[ flux.ghg$kappa.hist > max(kappa.hist)] <- max(kappa.hist) # higher end of simulation | |
| # flux.ghg$kappa.hist[ flux.ghg$kappa.hist < min(kappa.hist)] <- min(kappa.hist) # lower end of simulation | |
| # hist(flux.ghg$kappa.hist,40) | |
| lysi.gfg <- as.matrix(cbind(flux.ghg$flux.hist,flux.ghg$kappa.hist)) | |
| # plot(lysi.gfg[,1],lysi.gfg[,2],xlab="linear f0",ylab="kappa") | |
| # plot(log(lysi.gfg[,2]),lysi.gfg[,1],col="blue",ylab="estimated flux",xlab="kappa",pch=20,cex=0.5) | |
| # plot(log(flux.ghg$HMR.kappa),flux.ghg$rsq.rh,col="blue",ylab="estimated flux",xlab="kappa",pch=20,cex=0.5) | |
| # points(flux.ghg$HMR.f0,lysi.gfg[,2],col="blue") | |
| # | |
| # # lysi.gfg[order(lysi.gfg[,1]),] | |
| x <- abs(lysi.gfg[,1]) # vector with all fluxes | |
| y <- lysi.gfg[,2] # vector with all coresponding kappa | |
| # my_hist3d <- function(x, y, freq=FALSE, nclass="auto") { | |
| n<-length(x) | |
| kappa[1] <- 0.000001 | |
| # if (nclass == "auto") { nclass<-ceiling(sqrt(nclass.Sturges(x))) } | |
| nclass <- length(kappa) | |
| # breaks.x <- seq(min(x),max(x),length=(nclass+1)) | |
| # breaks.x <- seq(0,fmax,length.out=26)# c(f0,104)-2 | |
| breaks.x <- c(f0*0.95,fmax*1.05) # c(f0,104)-2 | |
| # breaks.y <- seq(min(y),max(y),length=(nclass+1)) | |
| breaks.y <- log(c(kappa*0.95,max(kappa)*1.05)) # c(0.0001,exp(seq(-11,-3.2,0.3))*3600) ; kappa = c(0.0001/36,exp(seq(-11,-3.5,0.32))) | |
| h <- NULL | |
| for (i in 1:nclass) | |
| for (j in 1:nclass) | |
| h <- c(h, sum(x <= breaks.x[j+1] & x >= breaks.x[j] & log(y) <= breaks.y[i+1] & log(y) >= breaks.y[i] ) ) | |
| # if (freq) h <- h / n | |
| # xx <-round(mean(breaks.x[1:2])+(0:(nclass-1))*diff(breaks.x[1:2]), 1) | |
| xx <-f0 | |
| xx[1] <- 0.000001 | |
| # yy <- as.factor(round(mean(breaks.y[1:2])+(0:(nclass-1))*diff(breaks.y[1:2]), 1)) | |
| # yy <- (length(breaks.y)-1):1 # as.factor(round(mean(breaks.y[1:2])+(0:(nclass-1))*diff(breaks.y[1:2]), 1)) | |
| yy <- 1:(length(breaks.y)-1) # as.factor(round(mean(breaks.y[1:2])+(0:(nclass-1))*diff(breaks.y[1:2]), 1)) | |
| res <- cbind(expand.grid(xx,yy), h) | |
| # kappa[1] <- 0.04/3600 | |
| colnames(res) <- c("f0","kappa","Frequency") | |
| res.matrix <- t(data.matrix(dcast(res,kappa~f0))[,-1]) | |
| res.matrix <- res.matrix[,-25] | |
| colnames(res.matrix) <- as.character(signif(kappa[-25],2)) | |
| # hist.results <- mean(flux.ghg[[fcs.method]],na.rm = T) # calculate the mean flux of the dataset | |
| hist.results <- geo_mean(flux.ghg[[fcs.method]]) # calculate the mean flux of the dataset | |
| hist.results[2] <- table(flux.ghg[[fcs.method]] == flux.ghg$HMR.f0 & !is.na(flux.ghg$HMR.kappa) )[2] # calculate the numer of HMR fluxes used | |
| hist.results[3] <- hist.results[1] / geo_mean(flux.ghg$linear.f0) # calcualte change from linear | |
| # fcs.parameter = c("IQR","bias","MSE") | |
| # for (i in 1:3) { | |
| ## calculate mean IQ90 | |
| iqr <- cbind(expand.grid(xx,kappa), simus.sum[simus.sum$method==fcs.method,IQR]) | |
| colnames(iqr) <- c("f0","kappa","IQR") | |
| iqr.matrix <- data.matrix(dcast(iqr,kappa~f0))[,-1] | |
| iqr.matrix <- iqr.matrix[,-25] | |
| iqr.matrix[,24] <- simus.sum$IQR[simus.sum$method=="linear.f0" & simus.sum$kappa > 2e-02] | |
| # iqr.matrix[14,18] <- 0.005 | |
| colnames(iqr.matrix) <- as.character(signif(kappa[-25],2)) | |
| ### frequency 2d histogramm ### | |
| # res.matrix[1,24] <- 400 | |
| tiff(filename = paste("graphs/Histplot-IQR_",fcs.method,".tif",sep=""), | |
| width = 1080, height = 800, units = "px", pointsize = 18) | |
| hist3D(z = res.matrix, x=1:25, y=log(kappa[-1]), border = "black", xlab = "log(flux)", ylab="log(kappa)", zlab="flux count", | |
| plot=T,add=F,colkey=T, zlim=c(0,160),theta=305, colvar= log(iqr.matrix), phi=45, resfac = 3, # zlim=c(0,450), log(SE.tot), # | |
| lighting=F,clab="log(IQ90)", clim = c(-7.5,4.5)) # ,col = ramp.col(c("darkgreen","yellow2","red"),alpha = 1)) # | |
| dev.off() | |
| hist.results[4] <- mean.param <- sum(res.matrix * iqr.matrix) / sum(res.matrix) | |
| # cbind(hist.results, ) | |
| # } | |
| ## calculate bias | |
| iqr <- cbind(expand.grid(xx,kappa), simus.sum[simus.sum$method==fcs.method,bias]) | |
| colnames(iqr) <- c("f0","kappa","bias") | |
| iqr.matrix <- data.matrix(dcast(iqr,kappa~f0))[,-1] | |
| iqr.matrix <- iqr.matrix[,-25] | |
| iqr.matrix[,24] <- simus.sum$bias[simus.sum$method=="linear.f0" & simus.sum$kappa > 2e-02] | |
| # iqr.matrix[14,18] <- 0.005 | |
| colnames(iqr.matrix) <- as.character(signif(kappa[-25],2)) | |
| ### frequency 2d histogramm ### | |
| # res.matrix[1,24] <- 400 | |
| tiff(filename = paste("graphs/Histplot-bias_",fcs.method,".tif",sep=""), | |
| width = 1080, height = 800, units = "px", pointsize = 18) | |
| hist3D(z = res.matrix, x=1:25, y=log(kappa[-1]), border = "black", xlab = "log(flux)", ylab="log(kappa)", zlab="flux count", | |
| plot=T,add=F,colkey=T, zlim=c(0,160),theta=325, colvar=iqr.matrix , phi=38, resfac = 3, # zlim=c(0,450), log(SE.tot), log(iqr.matrix+5) # | |
| lighting=F,clab="bias" ,col = ramp.col(c("red","orange","green","blue"),alpha = 1), clim = c(-5,0.5)) # ) # , clim = c(-3,2) | |
| dev.off() | |
| hist.results[5] <- mean.param <- sum(res.matrix * iqr.matrix) / sum(res.matrix) | |
| # cbind(hist.results, ) | |
| # } | |
| ## calcualte MSE | |
| iqr <- cbind(expand.grid(xx,kappa), simus.sum[simus.sum$method==fcs.method,MSE]) | |
| colnames(iqr) <- c("f0","kappa","MSE") | |
| iqr.matrix <- data.matrix(dcast(iqr,kappa~f0))[,-1] | |
| iqr.matrix <- iqr.matrix[,-25] | |
| iqr.matrix[,24] <- simus.sum$MSE[simus.sum$method=="kappa.max" & simus.sum$kappa <= 1.2e-05] | |
| # iqr.matrix[14,18] <- 0.005 | |
| colnames(iqr.matrix) <- as.character(signif(kappa[-25],2)) | |
| ### frequency 2d histogramm ### | |
| # res.matrix[1,24] <- 400 | |
| tiff(filename = paste("graphs/Histplot-MSE_",fcs.method,".tif",sep=""), | |
| width = 1080, height = 800, units = "px", pointsize = 18) | |
| hist3D(z = res.matrix, x=1:25, y=log(kappa[-1]), border = "black", xlab = "log(flux)", ylab="log(kappa)", zlab="flux count", | |
| plot=T,add=F,colkey=T, zlim=c(0,160),theta=325, colvar= log(iqr.matrix), phi=38, resfac = 3, # zlim=c(0,450), log(SE.tot), # | |
| lighting=F,clab="log(MSE)", clim = c(-16.3,3.2)) # ,col = ramp.col(c("darkgreen","yellow2","red"),alpha = 1)) # | |
| dev.off() | |
| hist.results[6] <- mean.param <- sum(res.matrix * iqr.matrix) / sum(res.matrix) | |
| # cbind(hist.results, ) | |
| # } | |
| hist.results[7] <- fcs.method | |
| write.table(t(hist.results), "hist_results.csv", sep = ";", row.names = F, col.names = F, append = T) | |
| } | |
| cols <- c("linear.f0","robust.linear.f0","HMR.f0","AIC","g.fact","kappa.max","RF2011","RF2017") | |
| # for (j in 1:8) hist.eval(flux.ghg, fcs.method = cols[j]) | |
| # hist.eval(flux.ghg, fcs.method ="robust.linear.f0") | |
| lapply(seq_along(cols), function(i) hist.eval(data = flux.ghg, cols[i])) | |
| # purrr::imap(cols, ~ hist.eval(data = flux.ghg, fcs.method = .x)) | |
| # | |
| library(plotly) | |
| p <- plot_ly(x = flux.ghg$flux.hist, y = log(flux.ghg$kappa.hist)) | |
| pp <- subplot( | |
| p %>% add_markers(alpha = 0.1), | |
| p %>% add_histogram2d() | |
| ) | |
| # colkey(clim = range(log(CI.tot2+0.0001)), clab = "log(90%IQR)", add = TRUE, dist = -0.05, length = 0.7, width = 0.7) | |
| # text3D(2, -5, 0, labels = expression(theta[1]), add = TRUE, adj = 1) | |
| # plotrgl(lighting = TRUE, smooth = TRUE) | |
| # | |
| # plot.gasfluxes.hist <- function(x, type = "hist", stat = "n") { | |
| # | |
| # if (type == "hist") { | |
| # | |
| # p <- ggplot(x, aes_(x = ~f0, y = ~kappa, z = as.name(stat))) + | |
| # stat_contour(aes(fill = ..level..), geom = "polygon", bins = 500) + | |
| # facet_wrap(~ method) + | |
| # scale_fill_distiller(name = sprintf("flux estimate\n(%s)", stat), palette = "RdYlGn") + # BrBG "YlGnBu" RdYlGn" | |
| # scale_y_continuous(trans = "log10", labels = function(b) signif(b, 3)) + | |
| # theme_bw() | |
| # return(print(p)) | |
| # | |
| # } | |
| # | |
| # } #more plot types can be implemented |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment