fix emiss*Rln error in cal_Rn

.gitignore

@@ -15,3 +15,4 @@ inst/doc
 
 Notes.md
 *.nc
+*.csv

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: PMLV2
 Title: Penman-Monteith-Leuning Evapotranspiration model V2
-Version: 0.1.0.9000
+Version: 0.1.0
 Authors@R: c(
   person("Dongdong", "Kong", role = c("aut", "cre"), 
          email = "[email protected]", comment = c(ORCID = "0000-0003-1836-8172")))
@@ -19,16 +19,17 @@ Imports:
     dplyr, lubridate, 
     data.table, 
     rtrend, terra, 
-    hydroTools, rtop,
+    hydroTools, 
     nctools,
     crayon
 Remotes:
     rpkgs/nctools, 
     CUG-hydro/hydroTools
 Suggests: 
     ncdf4,
-    knitr,
     pracma,
+    rtop,
+    knitr,
     rmarkdown,
     testthat (>= 3.0.0)
 Config/testthat/edition: 3

NAMESPACE

@@ -36,7 +36,6 @@ importFrom(hydroTools,cal_slope)
 importFrom(lubridate,yday)
 importFrom(lubridate,ymd)
 importFrom(purrr,map)
-importFrom(rtop,sceua)
 importFrom(stats,setNames)
 importFrom(terra,as.array)
 importFrom(terra,ext)

R/ET0_helper.R

@@ -1,7 +1,7 @@
 #' ET helper
 #' @name ET_helper
-#' 
-#' @description 
+#'
+#' @description
 #' - `get_lambda`: latent heat of vaporization (~2500 J g-1)
 #' - `get_Rn`: net solar radiation (W m-2)
 #' - `get_VPD`: vapor pressure deficit (kPa)
@@ -22,17 +22,17 @@ get_lambda <- function(Tavg) {
 #' @param Rln longwave inward solar radiation, W m-2
 #' @param albedo shortwave albedo, -
 #' @param emiss longwave emissivity, -
-#' 
+#'
 #' @rdname ET_helper
 #' @export
 get_Rn <- function(Rs, Rln, Tavg, albedo, emiss) {
   Stefan <- 4.903e-9 # Stefan-Boltzmann constant [MJ K-4 m-2 day-1],
 
   Rns <- (1 - albedo) * Rs
-  RLout <- emiss * Stefan * (Tavg + 273.15)^4
-  RLout <- RLout / 0.0864 # convert from MJ m-2 d-1 to W/m2
+  RLout <- Stefan * (Tavg + 273.15)^4 / 0.0864 # convert from MJ m-2 d-1 to W/m2
 
-  Rnl <- Rln - RLout
+  ## fix bug at emiss*Rln, 20230922
+  Rnl <- emiss*(Rln - RLout)
   Rn <- pmax(Rns + Rnl, 0) # Rns+Rnl, includenan
   Rn
 }

R/model_calib.R

@@ -1,12 +1,10 @@
 #' Calibration for `PML` model
 #'
 #' @inheritParams PML
-#' @inheritParams rtop::sceua
 #'
 #' @param IGBPcode (optional) [IGBP_006] code. If `hc_raw` specified in `par`, `IGBPcode` can be ignored.
 #' Otherwise, `IGBPcode` should be provided.
 #' 
-#' @importFrom rtop sceua
 #' @export
 PML_calib <- function(data, IGBPcode = 1, of_gof = NSE, ..., verbose = TRUE, maxn = 1e3) {
   # c("Alpha", "Thelta", "m", "Am_25", "D0", "kQ", "kA", "S_sls", "fER0", "VPDmin", "VPDmax")
@@ -26,7 +24,7 @@ PML_calib <- function(data, IGBPcode = 1, of_gof = NSE, ..., verbose = TRUE, max
   #     data = data, IGBPcode = IGBPcode, of_gof = of_gof)$optim$bestmem
   # par = RMPSH::RMPSH_opt(par0, PML_goal, lb, ub,
   #     data = data, IGBPcode = IGBPcode, of_gof = of_gof)
-  # par = sceua(PML_goal, par0, lb, ub, iprint = ifelse(verbose, 0, -1), maxn = maxn,
+  # par = rtop::sceua(PML_goal, par0, lb, ub, iprint = ifelse(verbose, 0, -1), maxn = maxn,
   #                 data = data, IGBPcode = IGBPcode, of_gof = of_gof)$par
   par <- set_names(par, parnames)
   par["hc"] <- options_param$hc_raw[IGBPcode]
@@ -123,6 +121,7 @@ cal_gof <- function(df, by = "IGBP") {
       dplyr::select(IGBP, R2, KGE, NSE, R, RMSE, Bias_perc, n_sim)
     ans
   }
+
   GPP <- df[, .(yobs = GPPobs, ysim = GPP, IGBP, site)] %>% .cal_gof()
   ET <- df[, .(yobs = ETobs, ysim = ET, IGBP, site)] %>% .cal_gof()
   info <- listk(ET, GPP) %>%

R/tidy_forcing_flux.R

@@ -45,7 +45,7 @@ tidy_forcing_flux <- function(df) # , par
 
   # 2. add VPD and Rn
   Rn <- get_Rn(df$Rs, df$Rl_in, Tavg, df$Albedo, df$Emiss)
-
+  
   # 2. intermediate variables
   rou_a <- 3.846 * 10^3 * Pa / (Tavg + 273.15) # kg m-3
   gama <- Cp * Pa / (0.622 * lambda) # kpa deg-1
@@ -55,7 +55,7 @@ tidy_forcing_flux <- function(df) # , par
   ## Equilibrium evaporation
   Eeq <- epsilon / (epsilon + 1) * Rn / lambda * 86400 * 10^-3 # mm
   Eeq <- pmax(0.0001, Eeq)
-
+  
   # Evaporation from soil surface
   if (is.null(df$f_value_soil)) {
     kA <- 0.9 # insensitive parameter for PMLv2

scripts/Fluxnet_calib_params/s1_calibrate_model_params.qmd

@@ -0,0 +1,90 @@
+```{r}
+# library(phenofit)
+library(hydroTools)
+library(Ipaper)
+library(data.table)
+library(dplyr)
+library(rfluxnet)
+# library(ggrepel)
+# library(ggpmisc)
+```
+
+## 1. 准备数据
+
+```{r}
+df_raw <- fread("data-raw/INPUTS/PMLv2_training_forcing_flux_v20200828 (80%)_102sp.csv")
+df_raw$date %<>% as.Date()
+
+# % considering Uz = 15 for fluxsites, and convert to real U2
+df <- df_raw[, .(site, IGBPcode, IGBPname, date, GPP_NT, GPP_DT, LE, LE_CORR,
+  LAI = LAI_sgfitw, LAI_sgfitw, LAI_raw,
+  dhour_norm,
+  Emiss,
+  CO2, Pa, Prcp, Rl_in, Rs, VPD, Tavg, U2 = cal_U2(U2, z.wind = 15), Albedo
+)]
+
+df %<>% tidy_forcing_flux()
+
+## 选择要率定参数的植被类型
+IGBPs <- c("CRO", "DBF", "EBF", "ENF", "GRA", "MF", "OSH", "SAV", "WET", "WSA")
+# IGBPs = c("EBF", "SAV", "WET")#[2]
+# IGBPs = "ENF"
+IGBPinfo <- df[, .(IGBPcode, IGBPname)] %>%
+  unique() %>%
+  .[match(IGBPs, IGBPname), ]
+n <- nrow(IGBPinfo)
+
+df_tmp <- df[IGBPname %in% IGBPs, .(site, IGBP = IGBPname, date, LAI_raw, LAI, LAI_sgfitw)]
+# pb <- progress::progress_bar$new(
+#     format = "[:bar] :percent eta: :eta, in :elapsed", total = n)
+```
+
+## 2. 参数率定
+
+```{r}
+set.seed(1) # 确保每次参数优选结果一致
+
+lst <- foreach(igbp_code = with(IGBPinfo, set_names(IGBPcode, IGBPname)), i = icount()) %do% {
+  runningId(i, prefix = IGBPinfo$IGBPname[i])
+
+  d_obs <- df[IGBPcode == igbp_code]
+  r <- PML_calib(d_obs, igbp_code, verbose = TRUE, maxn = 1e4)
+}
+
+res <- purrr::transpose(lst)[1:2] %>%
+  map(~ melt_list(map(.x, as.data.table), "IGBP"))
+
+res$data[, GOF(ETobs, ET)]
+
+res$data[, GOF(GPPobs, GPP)]
+```
+
+
+## 3. 拟合优度分析
+
+```{r}
+# res$gof %<>% .[match(IGBPs, IGBP), ]
+# # write_list2xlsx(res, "PMLV2_v1.xlsx")
+res_d8 <- merge(res$data, st_flux166 %>% rename(IGBP_org = IGBP),
+  all.x = TRUE, by = "site"
+) %>%
+  mutate(year = year(date) - 1 * (month(date) < 7) * (lat < 0)) # 考虑南半球的年份
+
+res_season <- res_d8 %>%
+  add_dn() %>%
+  .[, lapply(.SD, mean, na.rm = TRUE), .(IGBP, site, d8),
+    .SDcols = c("ET", "ETobs", "GPP", "GPPobs")
+  ]
+
+cat(green("[ALL]:"))
+gof <- list(d8 = res_d8, season = res_season) %>% map(cal_gof)
+gof
+# PML_gof(res_d8, of_gof = KGE) %>% {cat(str(.))}
+
+# write_list2xlsx(gof, "WHIT_Bisquare")
+# write_list2xlsx(gof, "SG_Bisquare_v0")
+
+# raw : 0.674, 0.602
+# SG  : 0.692, 0.75
+# WHIT: 0.682, 0.689
+```