add C3C4 task

Author: kongdd

Diff for: examples/data_PML_2019.jl

@@ -2,19 +2,6 @@ using Ipaper, HydroTools, PML
 using DataFrames, RTableTools
 using Test, BenchmarkTools
 
-function replace_miss!(df::AbstractDataFrame)
-  # colnames = names(df)
-  # num_cols = [name for name in colnames if getDataType(df[!, name]) <: Number]
-  for col in names(df)
-    x = df[!, col]
-    type = getDataType(x)
-    if type <: AbstractFloat
-      df[!, col] = drop_missing(x)
-    end
-  end
-  df
-end
-
 small = false
 if small
   f = "D:/GitHub/PML/PMLV2_Kong2019.m/data/INPUT_dt_v2019.csv"

Diff for: examples/ex2_NorthChina_C3C4_calib.jl

@@ -0,0 +1,38 @@
+using PML, Ipaper, Test, RTableTools
+
+include("main_pkgs.jl")
+data = fread("./Forcing_C3C4_sp2.csv")
+replace_miss!(data)
+
+df = data[data.site.=="固城", :]
+# df = data[data.site.=="禹城", :]
+# df.ET_obs = df.ET_obs ./ 0.9
+
+par = par0
+r = PMLV2(df; par)
+
+theta, goal, flag = model_calib(df, par0)
+df_out = PMLV2_sites(df; par=theta2par(theta))
+# df_out[1:10, :]
+sites = unique(df.site)
+
+
+begin
+  RES = []
+  for site in sites
+    inds = df.site .== site
+    d_obs = df[inds, :]
+    d_sim = df_out[inds, :]
+
+    gof = [
+      (; var="ET", GOF(d_obs.ET_obs, d_sim.ET)...),
+      (; var="GPP", GOF(d_obs.GPP_obs, d_sim.GPP)...)] |> DataFrame
+    gof = gof[:, [:var, :NSE, :KGE, :R2, :bias_perc, :n_valid]]
+    push!(RES, gof)
+  end
+  RES
+end
+
+## PML-ET
+# 固城: 0.46
+# 禹城: 0.63

Diff for: examples/main_pkgs.jl

@@ -0,0 +1,14 @@
+using DataFrames
+
+function replace_miss!(df::AbstractDataFrame)
+  # colnames = names(df)
+  # num_cols = [name for name in colnames if getDataType(df[!, name]) <: Number]
+  for col in names(df)
+    x = df[!, col]
+    type = getDataType(x)
+    if type <: AbstractFloat
+      df[!, col] = drop_missing(x)
+    end
+  end
+  df
+end

Diff for: src/Params.jl

@@ -19,8 +19,8 @@ $(TYPEDFIELDS)
   "initial slope of the CO2 response curve to assimilation rate, (i.e., carboxylation efficiency; `μmol m⁻² s⁻¹ [μmol m⁻² s⁻¹]⁻¹`)"
   Thelta::FT = 0.04  | (0.01, 0.07)
 
-  "stomatal conductance coefficient"
-  m     ::FT = 10.00 | (2.00, 100.00)
+  "stomatal conductance coefficient" # 气孔导度斜率参数
+  g1     ::FT = 10.00 | (2.00, 100.00)
 
   "carbon saturated rate of photosynthesis at 25 °C, `μmol m⁻² s⁻¹`"
   Am_25 ::FT = 50.00 | (5.00, 120.00)

Diff for: src/Utilize/model_calib.jl

@@ -15,11 +15,9 @@ function model_goal(df, theta; IGBPcode=nothing, of_gof=:NSE, verbose=false)
     indexes = [:KGE, :NSE, :R2, :RMSE, :bias, :bias_perc]
     info_GPP = info_GPP[indexes] |> round2
     info_ET = info_ET[indexes] |> round2
-    # @show info_GPP
-    # @show info_ET
   end
-
-  goal = (info_ET[of_gof] + info_GPP[of_gof]) / 2
+  gof = [info_ET[of_gof], info_GPP[of_gof]]
+  goal = weighted_mean(gof, [1, 0.5])
   goal
 end
 

Diff for: src/main_Ipaper.jl

@@ -61,3 +61,5 @@ end
 # export struct2vec, struct2tuple
 export round2;
 export movmean2, nanmean2, getDataType, replace_miss
+
+weighted_mean(x::AbstractVector, w::AbstractVector) = sum(x .* w) / sum(w)

Diff for: src/photosynthesis.jl

@@ -31,7 +31,7 @@ function photosynthesis(Tavg::T, Rs::T, VPD::T, LAI::T,
   GPP = Ag * 86400 / 10^6 * 12 # [umol m-2 s-1] to [g C m-2 d-1]
 
   f_VPD_gc = 1.0 / (1.0 + VPD / par.D0) # Leuning f_vpd
-  Gc = par.m * Ag / Ca * f_VPD_gc # canopy conductance for carbon
+  Gc = par.g1 * Ag / Ca * f_VPD_gc # canopy conductance for carbon
 
   ## Convert from mol m-2 s-1 to m s-1
   Gc = Gc * 1e-2 / (0.446 * (273 / (273 + Tavg)) * (Pa / 101.3)) # Gc = Gc * mol2m(Tavg, Pa)