second test whittaker lambda formula

Author: kongdd

DESCRIPTION

@@ -8,7 +8,7 @@ Description: A remote sensing vegetation index phenology extraction package.
 License: MIT | file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 6.0.1.9000
+RoxygenNote: 6.0.1
 LinkingTo: Rcpp, RcppArmadillo
 Imports: Rcpp, RcppArmadillo, 
     tibble, dplyr, purrr, stringr, tidyr, ggplot2, lubridate, 

NEWS.md

@@ -0,0 +1,6 @@
+# phenofit 0.1.2
+
+* Added a `NEWS.md` file to track changes to the package.
+
+
+

R/tools.R

@@ -199,22 +199,39 @@ cv_coef <- function(x, w){
 #' @param Y_sim Numeric vector, corresponding simulated values
 #' @param w Numeric vector, weights of every points
 #' @param include.cv If true, cv will be returned. 
+#' 
+#' @return 
+#' \itemize{
+#' \item \code{RMSE} root mean square error
+#' \item \code{NSE} NASH coefficient
+#' \item \code{R2} correlation of determination
+#' \item \code{Bias} bias
+#' \item \code{Bias_perc} bias percentage 
+#' \item \code{MAE} mean absolute error 
+#' \item \code{pvalue} pvalue of \code{R}
+#' \item \code{n_sim} number of valid obs.
+#' \item \code{R} pearson correlation
+#' }
 #' @export
 GOF <- function(Y_obs, Y_sim, w, include.cv = FALSE){
     if (missing(w)) w <- rep(1, length(Y_obs))
 
-    # remove NA values in Y_sim, Y_obs and w
-    I <- which(!(is.na(Y_sim) | is.na(Y_obs) | is.na(w)))
+    # remove NA and Inf values in Y_sim, Y_obs and w
+    valid <- function(x) !is.na(x) & is.finite(x) 
+    
+    I <- which(valid(Y_sim) & valid(Y_obs) & valid(w))
     # n_obs <- length(Y_obs)
     n_sim <- length(I)
 
     Y_sim <- Y_sim[I]
     Y_obs <- Y_obs[I]
+    w     <- w[I]
 
     if (include.cv) CV <- cv_coef(Y_obs, w)
     if (is_empty(Y_obs)){
-        out <- c(Bias = NA, MAE = NA,RMSE = NA, NSE = NA, R2 = NA,
-             pvalue = NA, n_sim = NA, R = NA)
+        out <- c(RMSE = RMSE, NSE = NSE, R2 = R2, MAE = MAE, 
+            Bias = Bias, Bias_perc = Bias_perc, 
+            R = NA, pvalue = NA, n_sim = NA)
         if (include.cv) out <- c(out, CV)
         return(out) #R = R,
     }
@@ -228,10 +245,12 @@ GOF <- function(Y_obs, Y_sim, w, include.cv = FALSE){
 
     SSR    <- sum( (Y_sim - y_mean)^2 * w) 
     SST    <- sum( (Y_obs - y_mean)^2 * w)
-    R2     <- SSR / SST
+    # R2     <- SSR / SST
+    R2     <- summary(lm(Y_sim ~ Y_obs))$r.squared
 
     RE     <- Y_sim - Y_obs
     Bias   <- sum ( w*RE)     /sum(w)                     # bias
+    Bias_perc <- Bias/y_mean                              # bias percentage        
     MAE    <- sum ( w*abs(RE))/sum(w)                     # mean absolute error
     RMSE   <- sqrt( sum(w*(RE)^2)/sum(w) )                # root mean sqrt error
 
@@ -251,13 +270,15 @@ GOF <- function(Y_obs, Y_sim, w, include.cv = FALSE){
         message(e$message)
     })
 
-    out <- c(Bias = Bias, MAE = MAE,RMSE = RMSE, NSE = NSE, R2 = R2,
-             pvalue = pvalue, n_sim = n_sim, R = R)
+    out <- c(RMSE = RMSE, NSE = NSE, R2 = R2, MAE = MAE, 
+             Bias = Bias, Bias_perc = Bias_perc, 
+             R = R, pvalue = pvalue, n_sim = n_sim)
     if (include.cv) out <- c(out, CV)
     return(out)
 }
 
 
+
 #' Determinated correlation critical value
 #' 
 #' @param n length of observation.

man/GOF.Rd

@@ -19,6 +19,7 @@ if (!file.exists(file)){
 }else{
     load(file)
 }
+
 # load("D:/Documents/GoogleDrive/phenofit.rda")
 
 methods <- c('AG', 'BECK', 'ELMORE', 'ZHANG', 'whit_R', 'whit_gee')[-5]
@@ -29,35 +30,58 @@ prefix <- c("phenoflux", "phenocam")[i]
 df <- if(i == 1) df_flux else df_cam
 st <- if(i == 1) st_flux else st_cam
 
-df      <- df[iters == "iter2"]
 st$IGBPname %<>% factor(IGBPnames_006)
+st <- st[order(IGBPname,site), ] # reorder according to IGBP
+st[site %in% sites_sel, titlestr := sprintf("(%s) %s, %s", letters[1:.N],site, IGBPname)]
 
 # make sure different curve fitting methods have the same length fitting
-formula <- if(i == 1) formula(site+date+t+y+GPP_NT+GPP_DT+SummaryQA~meth) else
-    formula(site+date+t+y+gcc+vci+SummaryQA~meth)
-
-over_perform(df, formula, prefix)
+formula <- if(i == 1) formula(site+date+t+y+GPP_NT+GPP_DT+SummaryQA+iters~meth) else
+    formula(site+date+t+y+gcc+vci+SummaryQA+iters~meth)
+IGBP.all = F
+if (i == 1){
+    over_perform(df[iters == "iter2"], formula, prefix, IGBP.all = IGBP.all)
+} else{
+    over_perform(df[iters == "iter2"], formula, prefix, ylim2 = c(64, 100), IGBP.all = IGBP.all)
+}
 
 # site figure data input
-
 df_trim <- dcast(df, formula, value.var = "value", fun.aggregate = mean) # %>% na.omit()
 df_trim$SummaryQA %<>% factor(qc_levels)
 # df_trim <- melt(df_trim, measure.vars = methods, variable.name = "meth")
 
-sites <- unique(df$site)
+sites    <- st$site
+# sites <- unique(df$site)
 sitename <- sites[100]
 
-vars <- c("get_range", "save_pdf", "lgd_vci", "lgd_gpp",
-          "qc_levels", "qc_colors", "qc_shapes", "methods")
-cl <- cluster_Init(pkgs = c("data.table", "ggplot2", "magrittr"))
-clusterExport(cl, vars)
-res <- parLapplyLB(cl, sites, plot_whit,
-                   df_trim, st, prefix_fig = paste0("whit_", prefix))
-for (i in seq_along(sites)){
-    runningId(i)
-    sitename <- sites[i]
-    plot_whit(sitename, df_trim, st, prefix_fig = paste0("whit_", prefix))
+# vars <- c("get_range", "save_pdf", "lgd_vci", "lgd_gpp",
+#           "qc_levels", "qc_colors", "qc_shapes", "methods")
+# cl <- cluster_Init(pkgs = c("data.table", "ggplot2", "magrittr"))
+# clusterExport(cl, vars)
+# res <- parLapplyLB(cl, sites, plot_whit,
+#                    df_trim, st, prefix_fig = paste0("whit_", prefix))
+
+plot_whits <- function(sites, df_trim, file){
+    ps <- list()
+    for (i in seq_along(sites)){
+        runningId(i)
+        sitename <- sites[i]
+        # for single method, ggplot obj return
+        ps[[i]] <- plot_methods(sitename, df_trim, st, prefix_fig = paste0("whit_", prefix), methods = "whit_gee", show.legend = F)
+        # print(p)
+        # for all methods, grob obj return
+        # p <- plot_methods(sitename, df_trim, st, prefix_fig = paste0("whit_", prefix), methods)
+    }
+
+    ylab_r  <- expression("GPP ( gC "*mm^-1*d^-1*" )")
+    # ylab_r  <- "VCI"
+    FigsToPages(ps, lgd_gpp, ylab_r, file, width = 10)
 }
+plot_whits(sites, df_trim, "gee_whit_flux166.pdf") # all sites
+
+## 2. select representative points for fluxnet
+sites_sel <- c("RU-Fyo", "FR-Pue", "US-UMB", "BE-Vie", "US-KS2", "US-Whs",
+    "AU-How", "AU-Dry", "CH-Fru", "US-Los", "DE-Geb")
+plot_whits(sites_sel, df_trim, "gee_whit_flux11.pdf") # all sites
 
 # merge_pdf('../whit_phenoflux166.pdf', indir = "./")
 # merge_pdf('whit_phenoflux166.pdf', indir = "Figure/")
@@ -77,6 +101,5 @@ merge_pdf('whit_phenoflux166.pdf', indir = "Figure", "whit_phenoflux.*.pdf", del
 #                                   "snow/ice" = "#F8766D", "cloud" = "#C77CFF"), drop = F) +
 #     scale_shape_manual(values = c(19, 15, 4, 17), drop = F) +
 #     scale_y_continuous(lim = lim_raw)
-
 # p3 <- ggplot_dual_axis(p1, p2) #%>% as.ggplot()
 # p <- ggplot_dual_axis(p3, p0, add_yaxis_r = F)

test/07_whit/03_evaluate_Whittaker_GOF.R

@@ -26,7 +26,7 @@ lgd   <- phenofit:::make_legend(linename = c("iter1", "iter2"),
                                 linecolor = c("blue", "red"))
 
 file <- "whit_GEE.pdf"
-Cairo::CairoPDF(file, 10, 4)
+# Cairo::CairoPDF(file, 10, 4)
 # par(mfrow = c(4, 1), mar = c(1, 2, 3, 1), mgp = c(1.5, 0.6, 0))
 ## sometimes sample and reduceRegions result is different
 ps  <- list()
@@ -64,30 +64,7 @@ for (i in seq_along(sites)){
     # grid.draw(p2)
 }
 
-dev.off()
-file.show(file)
+# dev.off()
+# file.show(file)
 
 file = "Fig3_whit_point_example.pdf"
-Cairo::CairoPDF(file, 10, nrow*1.6)
-
-nrow  <- 6
-npage <- ceiling(length(ps)/nrow)
-for (i in 1:npage){
-    runningId(i)
-
-    I_beg <- (i - 1) * nrow + 1
-    I_end <- min(i*nrow, length(ps))
-
-    x  <- c(ps[I_beg:I_end], list(lgd))
-    nx <- length(x)
-    x[[nx - 1]] <- x[[nx - 1]] +
-        theme(axis.title.x = element_text(),
-              axis.text.x = element_text())
-    p2 <- gridExtra::arrangeGrob(grobs = x, nrow = nx, ncol = 1,
-                                 heights = c(rep(5, nx - 2), 5.5, 1),padding = unit(1, "line"),
-                                 left = textGrob("EVI", rot = 90, gp=gpar(fontsize=14)))
-    if (i != 1) grid.newpage();
-    grid.draw(p2)
-}
-dev.off()
-file.show(file)