PMLV2.jl passed test

Author: kongdd

Diff for: data/CRO/FLUXNET_CRO

@@ -1,4 +1,3 @@
-
 function calib_PML(data::AbstractDataFrame; of_gof=:KGE, maxn=2500)
   vars = ["IGBPname", "IGBPcode", "site", "date", "GPP_obs", "ET_obs",
     "Prcp", "Tavg", "U2", "Rn", "Rs", "VPD", "LAI", "Pa", "Ca"]
@@ -33,7 +32,7 @@ end
 # data.LAI .= data.LAI_whit
 
 # data = @rename(data, LAI = LAI_whit) # LAI_raw, LAI_whit, LAI_sgfitw
-r = calib_PML(data; of_gof=:KGE, maxn=500)
+r = calib_PML(data; of_gof=:KGE, maxn=5000)
 r.gof
 
 # out = fread("D:/GitHub/PML/PMLV2_Kong2019.m/OUTPUT/PMLv2_flux102_Cal_flux_v012.csv")

Diff for: scripts/calib_PML.jl

@@ -18,3 +18,28 @@ function aerodynamic_conductance(U2::T, hc::T; Zob=15.0) where {T<:Real}
   Ga
 end
 
+cal_rho_a(Tair, Pa) = 3.486 * Pa / 1.01(Tair + 273.15) # kg/m3
+# cal_rho_a(Tair, Pa) = 3.846 * Pa / (Tair + 273.15) # [kg/m3], error v2019
+
+# # rho_a: kg m-3
+# function cal_rho_a(Tair, q, Pa)
+#   3.486 * Pa / cal_TvK(Tair, q) # FAO56, Eq. 3-5
+# end
+
+# cal_rho_a(Tair, Pa) = 3.486 * Pa / cal_TvK(Tair)
+# # cal_rho_a(Tair, Pa) = 3.486 * Pa / cal_TvK(Tair)
+
+# # FAO56, Eq. 3-7
+# cal_TvK(Tair) = 1.01 * (Tair + 273)
+
+# # https://github.com/CUG-hydro/class2022_CUG_HydroMet/blob/master/ch02_大气的基本特征.md
+# # q ≈ ϵ*ea/Pa
+# # q = ϵ*ea/(Pa - (1 - ϵ)*ea)
+# function cal_TvK(Tair, q)
+#   # ea / Pa = q/(ϵ + (1 - ϵ) * q)
+#   (Tair + K0) * (1 + (1 - epsilon) * q / (ϵ + (1 - ϵ) * q))
+# end
+
+# # 这个是最精确的版本
+# # FAO56, Eq. 3-6
+# cal_TvK(Tair, ea, Pa) = (Tair + K0) * (1 + (1 - epsilon) * ea / Pa)

Diff for: src/ET_helper.jl

@@ -2,10 +2,11 @@ module PML
 
 export PMLV2_sites
 export PMLV2, T_adjust_Vm25, f_VPD_Zhang2019
+export GOF
 
 using DocStringExtensions
 using Parameters, DataFrames
-import HydroTools: cal_rho_a, cal_Uz, ET0_eq, 
+import HydroTools: cal_Uz, ET0_eq, 
   Cp, atm, 
   GOF, sceua
 

Diff for: src/PML.jl

@@ -2,7 +2,7 @@ using Statistics
 
 round2(x::NamedTuple, digits=3; kw...) = map(val -> round(val; digits), x)
 
-# main script of moving average
+
 function movmean2(y::AbstractVector{T}, win_left::Integer, win_right::Integer=win_right) where {T<:Real}
   n = length(y)
   z = zeros(Float64, n)
@@ -11,15 +11,14 @@ function movmean2(y::AbstractVector{T}, win_left::Integer, win_right::Integer=wi
     i_beg = max(i - win_left, 1)
     i_end = min(i + win_right, n)
 
-    n_i = 0        # number
-    sum = T(0.0)   # sum of values in window
-
+    count = 0    # number
+    ∑ = T(0.0)   # sum of values in window
     for j in i_beg:i_end
       isnan(y[j]) && continue # skip missing values
-      n_i += 1
-      sum += y[j]
+      count += 1
+      ∑ += y[j]
     end
-    z[i] = n_i > 0 ? sum / n_i : T(NaN)
+    z[i] = count > 0 ? ∑ / count : T(NaN)
   end
   z
 end

Diff for: src/main_Ipaper.jl

@@ -55,6 +55,7 @@ function PMLV2(Prcp::T, Tavg::T, Rs::T, Rn::T, VPD::T, U2::T, LAI::T,
   ρa = cal_rho_a(Tavg, Pa)
   # ρa = cal_rho_a(Tavg, q, Pa)
   LEca = (ρa * Cp * 1e6 * r.Ga * VPD / γ) / (ϵ + 1 + r.Ga / r.Gc_w) # W m-2, `Cp*1e6`: [J kg-1 °C-1]
+  # [kg m-3] [J kg-1 K-1] [m s-1] [kPa] / [kPa K-1] = [W m-2]
 
   r.Ecr = W2mm(LEcr; λ) # [W m-2] change to [mm d-1]
   r.Eca = W2mm(LEca; λ) # [W m-2] change to [mm d-1]
@@ -70,18 +71,6 @@ function PMLV2(Prcp::T, Tavg::T, Rs::T, Rn::T, VPD::T, U2::T, LAI::T,
 end
 
 
-"""
-# Arguments
-- `kw`: named keyword arguments
-  + `r`: `interm_PML`
-"""
-function PMLV2(d::AbstractDataFrame; par::Param_PMLV2=Param_PMLV2(), kw...)
-  PMLV2(d.Prcp, d.Tavg, d.Rs, d.Rn,
-    d.VPD, d.U2, d.LAI,
-    d.Pa, d.Ca; par, kw...) |> to_df
-end
-
-
 """
     PMLV2(Prcp, Tavg, Rs, Rn, VPD, U2, LAI, Pa, Ca; par=param0, frame=3)
 
@@ -119,3 +108,15 @@ function PMLV2(Prcp::V, Tavg::V, Rs::V, Rn::V,
   res.ET .= res.Ec .+ res.Ei .+ res.Es
   res
 end
+
+
+"""
+# Arguments
+- `kw`: named keyword arguments
+  + `r`: `interm_PML`
+"""
+function PMLV2(d::AbstractDataFrame; par::Param_PMLV2=Param_PMLV2(), kw...)
+  PMLV2(d.Prcp, d.Tavg, d.Rs, d.Rn,
+    d.VPD, d.U2, d.LAI,
+    d.Pa, d.Ca; par, kw...) |> to_df
+end

Diff for: src/model_PMLV2.jl

@@ -1,15 +1,49 @@
-using PML, Test
+using PML, Test, Ipaper
 
-@testset "ET PMLV2" begin
+function build_data()
+  param = fread("data/CRO/PARAM.csv")
+  df_out = fread("data/CRO/OUTPUT.csv")
+  df_out.Ec = df_out.Ecr + df_out.Eca
+
+  df = fread("data/CRO/INPUT.csv")
+  df.Ca = df.CO2
+
+  inds = [1:4; 10:11; 5:9]
+  _names = names(param)[inds]
+  _values = Matrix(param)[inds]
+  par = Param_PMLV2(_values..., 0.5)
+
+  l = (; df_out, df, par)
+  serialize("data/CRO/FLUXNET_CRO", l)
+end
+
+dir_proj = "$(@__DIR__)/.."
+df_out, df, par = deserialize("$dir_proj/data/CRO/FLUXNET_CRO")
+r = PMLV2_sites(df; par)
+
+
+@testset "PMLV2 scalar" begin
   Prcp = 2.0 # mm
   Tavg = 20.0
   Rs = 200.0
   Rn = 50.0
   VPD = 2.0
   U2 = 2.0
   LAI = 2.0
-  
-  # par = param0
-  # par = add(par, list(hc=2.0))
   @test_nowarn PMLV2(Prcp, Tavg, Rs, Rn, VPD, U2, LAI;)
 end
+
+@testset "CHECK PMLV2 RESULT" begin
+  @test GOF(df.Eeq, r.Eeq).MAE <= 0.01
+  @test GOF(df_out.Es_eq, r.Es_eq).MAE <= 0.002
+  @test GOF(df_out.Ga, r.Ga).MAE <= 1e-10
+  @test GOF(df_out.Gc, r.Gc_w).MAE <= 1e-3
+  @test GOF(df_out.Ei, r.Ei).MAE <= 1e-8 # Ei passed Test
+  @test GOF(df_out.Ec, r.Ec).MAE <= 0.002
+  @test GOF(df_out.Es, r.Es).MAE <= 0.01
+  @test GOF(df_out.Ecr, r.Ecr).MAE <= 0.002
+  @test GOF(df_out.Eca, r.Eca).MAE <= 0.002
+
+  @test GOF(df_out.ET_sim, r.ET).MAE <= 0.015
+  @test GOF(df_out.GPP_sim, r.GPP).MAE <= 1E-8
+end