minor update

Project.toml

@@ -10,6 +10,7 @@ LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"

TODO.md

@@ -2,8 +2,6 @@
 
 > 20230924
 
-- [x] 1. 检查计算结果是否正确
+- [ ] 1. 检查`PMLV2`敏感性参数
 
-- [ ] 2. 添加参数率定模块
-
-1. 把`Ca`添加到df
+- [ ] 2. 测试LAI除噪的效果

scripts/ET/s1_calibrate_model_params.qmd

@@ -46,10 +46,11 @@ vars = ["IGBPname", "IGBPcode", "site", "date", "GPP_obs", "ET_obs",
 
 IGBPs = unique_sort(data.IGBPname)
 IGBP = IGBPs[1]
-df = data[data.IGBP .== IGBP, vars]
+df = data[data.IGBP.==IGBP, vars]
 
 # @time theta, goal, flag = m_calib(df; IGBPcode, maxn=2500);
 @time theta, goal, flag = m_calib(df; IGBPcode=df.IGBPcode[1], maxn=2500);
+
 m_goal(df, theta; verbose=true)
 
 r = PMLV2_sites(df; par=theta2par(theta))
@@ -63,30 +64,35 @@ cbind(df[:, [:site, :date, :ET_obs, :GPP_obs]], r)
 
 ```{julia}
 res = []
+params = []
 
 for i in eachindex(IGBPs)
   IGBP = IGBPs[i]
   df = data[data.IGBP .== IGBP, vars]
-  printstyled("[i=$i, IGBP = $(IGBP)] -------------------------------\n", 
-    bold=true, color=:green, underline=false)
 
-  @time theta, goal, flag = m_calib(df; IGBPcode=df.IGBPcode[1], maxn=2500);
-  
-  printstyled("-----------------------------\n", color=:blue)
+  @time theta, goal, flag = m_calib(df; IGBPcode=df.IGBPcode[1], maxn=2500)
+
+  printstyled("[i=$i, IGBP = $(IGBP)] -------------------------------\n",
+    bold=true, color=:green, underline=false)
   m_goal(df, theta; verbose=true)
+  println("-----------------------------------------------------------\n")
 
+  par = theta2par(theta)
   ## 输出数据
-  r = PMLV2_sites(df; par=theta2par(theta))
+  r = PMLV2_sites(df; par)
   r = cbind(df[:, [:site, :date, :ET_obs, :GPP_obs]], r)
 
-  ## 需要记录结果，记录参数
   push!(res, r)
+  push!(params, par)
 end
 
 df_out = vcat(res...)
 
 GOF(df_out.ET_obs, df_out.ET)
 GOF(df_out.GPP_obs, df_out.GPP)
+
+mat_param = map(collect, params) |> x -> cat(x..., dims=2) |> transpose
+d_param = DataFrame(mat_param, parNames) |> d -> cbind(DataFrame(IGBP=IGBPs), d)
 ```
 
 

src/ET/calibrate.jl

@@ -59,9 +59,11 @@ end
 
 
 function m_goal(df, theta; IGBPcode=nothing, of_gof=:NSE, verbose=false)
+  # IGBPcode !== nothing && (par.hc = hc_raw[IGBPcode])
+  IGBPcode !== nothing && (theta[end] = hc_raw[IGBPcode]) # the last one is hc
   par = list(parNames, theta)
-  IGBPcode !== nothing && (par.hc = hc_raw[IGBPcode])
-
+  # @show theta
+  
   dobs = df[!, [:GPP_obs, :ET_obs]]
   dsim = PMLV2_sites(df; par)
 
@@ -86,7 +88,7 @@ end
 
 ## 最后一步，参数率定模块
 function m_calib(df::AbstractDataFrame; IGBPcode=nothing, maxn=2500, kw...)
-  theta, goal, flag = sceua(theta -> -m_goal(df, theta; kw...),
+  theta, goal, flag = sceua(theta -> -m_goal(df, theta; IGBPcode, kw...),
     theta0, parRanges[:, 1], parRanges[:, 2]; maxn, kw...)
   theta, goal, flag
 end

src/Optim/sceua.jl

@@ -12,6 +12,7 @@ x, feval, exitflag, output = sceua(fn, x0, bl, bu)
 ```
 """
 function sceua(fn::Function, x0::Vector, bl::Vector, bu::Vector;
+  verbose=false,
   maxn=1000, kstop=5, pcento=0.01, peps=0.0001, ngs=5, iseed=1, iniflg=1)
   ## This is the subroutine implementing the SCE algorithm
   # written by Q.Duan; 9/2004
@@ -230,12 +231,14 @@ function sceua(fn::Function, x0::Vector, bl::Vector, bu::Vector;
     # Check for convergency
     if icall >= maxn
       exitflag = 0
-      disp("\n*** Optimization search terminated because the limit ***")
-      println("On the maximum number of trials $(maxn) has been exceeded!")
+      if verbose
+        disp("\n*** Optimization search terminated because the limit ***")
+        println("On the maximum number of trials $(maxn) has been exceeded!")
+      end
     end
     if gnrng .< peps
       exitflag = 1
-      disp("The population has converged to a prespecified small parameter space")
+      verbose && disp("The population has converged to a prespecified small parameter space")
     end
     push!(criter, bestf)
     # criter = [criter bestf]'
@@ -245,16 +248,19 @@ function sceua(fn::Function, x0::Vector, bl::Vector, bu::Vector;
 
       if criter_change .< pcento
         exitflag = 1
-        println("The best point has improved in last $(num2str(kstop)) loops by less than the threshold $(num2str(pcento))")
-        println("Convergency has achieved based on objective function criteria!!!")
+        if verbose
+          println("The best point has improved in last $(num2str(kstop)) loops by less than the threshold $(num2str(pcento))")
+          println("Convergency has achieved based on objective function criteria!!!")
+        end
       end
     end
     # End of the Outer Loops
   end
 
-  @printf("Search was stopped at trial number: %d \n", icall)
-  println("Normalized geometric range = $(num2str(gnrng))")
-  println("The best point has improved in last $(num2str(kstop)) LOOPS BY $(num2str(criter_change))")
-
+  if verbose
+    @printf("Search was stopped at trial number: %d \n", icall)
+    println("Normalized geometric range = $(num2str(gnrng))")
+    println("The best point has improved in last $(num2str(kstop)) LOOPS BY $(num2str(criter_change))")
+end
   bestx, bestf, exitflag
 end

src/cal_radiation.jl

@@ -1,7 +1,7 @@
 export cal_Rn, cal_Rsi, cal_Rsi_toa, 
   cal_Rln, cal_Rln_out, cal_Rli, cal_Rln_yang2019
 
-
+# 辐射用的单位都是W
 """
     cal_Rn(Rs, Rln, Tavg, albedo, emiss)
 
@@ -30,9 +30,7 @@ function cal_Rn(Rs::T, Rln::T, Tavg::T, α::Float64, ϵ::Float64) where {T<:Real
   # max(Rn, T(0.0))
 end
 
-# convert from MJ m-2 d-1 to W/m2
-cal_Rln_out(T::FT, ϵ=0.96) where {FT<:Real} = Stefan * (T+K0)^4 / 0.0864 
-
+cal_Rln_out(T::FT, ϵ=1.0) where {FT<:Real} = ϵ * σ * (T + K0)^4
 
 """
   cal_Rsi_toa(lat, J)
@@ -54,14 +52,14 @@ function cal_Rsi_toa(lat=0, J::Integer=1)
   # 24 * 60 * 0.082 = 118.08
   lat = deg2rad(lat)
   Rsi_toa = 118.08 * dr / π * (ws * sin(lat) * sin(σ) + cos(lat) * cos(σ) * sin(ws)) # Allen, Eq. 21
-  max(Rsi_toa, 0)
+  max(MJ2W(Rsi_toa), 0)
 end
 
 
 """
   cal_Rsi(lat, J, ssd=nothing, cld=nothing, Z=0, a=0.25, b=0.5)
 
-Daily inward shortwave solar radiation at crop surface in MJ m-2 day-1 by
+Daily inward shortwave solar radiation at crop surface in W m-2 by
 providing sunshine duration (SSD) in hours or cloud cover in fraction.
 
 # Arguments
@@ -78,7 +76,7 @@ providing sunshine duration (SSD) in hours or cloud cover in fraction.
 - `b`: Coefficient of the Angstrom formula. Default is 0.50.
 
 # Returns
-- A tuple of solar radiation at crop surface in MJ m-2 day-1:
+- A tuple of solar radiation at crop surface in W m-2:
   - `Rsi`: Surface downward shortwave radiation.
   - `Rsi_o`: Clear-sky surface downward shortwave radiation.
   - `Rsi_toa`: Extraterrestrial radiation.
@@ -128,7 +126,6 @@ Net outgoing longwave radiation.
 """
 function cal_Rln(Tmax, Tmin, ea, cld)
   cld = isnan(cld) ? 1.0 : cld
-
   return (4.093e-9 * (((Tmax + 273.15)^4 + (Tmin + 273.15)^4) / 2)) *
          (0.34 - (0.14 * sqrt(ea))) *
          (1.35 * (1.0 - cld) - 0.35)

src/unit_convert.jl

@@ -19,11 +19,11 @@ R2Q(R::Real, area::Real; dt=24) = R * area / (dt * 3.6)
 Q2R(Q::Real, area::Real; dt=24) = Q * dt * 3.6 / area
 
 
-MJ2W(x::Real) = x / 86400 * 1e6
+MJ2W(x::Real) = x / 0.0864 # x / 86400 * 1e6
 
 MJ2mm(x::Real) = x / 2.45
 
-W2MJ(x::Real) = x / 1e6 * 86400
+W2MJ(x::Real) = x * 0.0864
 
 function W2mm(x::Real, Tair=0)
   lamada = 2500 - 2.2 * Tair