From 847103edc5455e927218baa15f24b895a8120634 Mon Sep 17 00:00:00 2001 From: Dongdong Kong Date: Sun, 22 Dec 2024 12:38:56 +0800 Subject: [PATCH] =?UTF-8?q?=E6=94=B9=E8=BF=9B=E5=8F=82=E6=95=B0=E4=BC=98?= =?UTF-8?q?=E5=8C=96=E6=96=B9=E6=A1=88?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/{Canopy => }/Canopy.jl | 0 src/{Optim => }/DataType.jl | 12 +++++ src/ModelCalib.jl | 95 +++++++++++++++++++++++++++++++++++++ src/Optim/ModelCalib.jl | 54 --------------------- src/Optim/model_gof.jl | 56 ---------------------- src/PML.jl | 18 ++++--- src/Params.jl | 25 +++++----- test/test-PMLV2.jl | 38 ++++++++------- 8 files changed, 152 insertions(+), 146 deletions(-) rename src/{Canopy => }/Canopy.jl (100%) rename src/{Optim => }/DataType.jl (85%) create mode 100644 src/ModelCalib.jl delete mode 100644 src/Optim/ModelCalib.jl delete mode 100644 src/Optim/model_gof.jl diff --git a/src/Canopy/Canopy.jl b/src/Canopy.jl similarity index 100% rename from src/Canopy/Canopy.jl rename to src/Canopy.jl diff --git a/src/Optim/DataType.jl b/src/DataType.jl similarity index 85% rename from src/Optim/DataType.jl rename to src/DataType.jl index 8ef5ff3..88d1523 100644 --- a/src/Optim/DataType.jl +++ b/src/DataType.jl @@ -75,6 +75,18 @@ function to_df(res::output_PML{T}) where {T<:Real} DataFrame(data, names) end +function map_df_tuple(fun::Function, lst::GroupedDataFrame{DataFrame}, args...; kw...) + n = length(lst) + _keys = keys(lst) + map(i -> begin + d = lst[i] + key = NamedTuple(_keys[i]) + r = fun(d, args...; kw...) + (; key..., r...) + end, 1:n) +end + export interm_PML, output_PML export to_mat; +export map_df_tuple diff --git a/src/ModelCalib.jl b/src/ModelCalib.jl new file mode 100644 index 0000000..fd73819 --- /dev/null +++ b/src/ModelCalib.jl @@ -0,0 +1,95 @@ +export ModelCalib, ModelCalib_IGBPs +export model_goal +export model_gof + +include("DataType.jl") +# include("model_gof.jl") + +function ModelCalib_IGBPs(data::AbstractDataFrame; + parNames=ParNames, + of_gof=:KGE, maxn=2500) + + vars = ["IGBPname", "IGBPcode", "site", "date", "GPP_obs", "ET_obs", + "Prcp", "Tavg", "U2", "Rn", "Rs", "VPD", "LAI", "Pa", "Ca"] + IGBPs = unique(data.IGBPname) |> sort + + @time params = par_map(IGBP -> begin + df = data[data.IGBP.==IGBP, vars] + IGBPcode = df.IGBPcode[1] + theta, _, _ = ModelCalib(df, par0, parNames; IGBPcode, of_gof, maxn, verbose=false) + theta + end, IGBPs; progress=false) + + # printstyled("[i=$i, IGBP = $IGBP] \n", bold=true, color=:green, underline=false) + res = map(i -> begin + IGBP = IGBPs[i] + theta = params[i] + df = data[data.IGBP.==IGBP, vars] + model_goal(df, theta, parNames; verbose=true) # print GOF + + par = theta2par(theta, parNames) + r = PMLV2_sites(df; par) + cbind(df[:, [:site, :date, :ET_obs, :GPP_obs]], r) + end, eachindex(params)) + df_out = melt_list(res; IGBP=IGBPs) + + gof = model_gof(df_out) + (; param=theta2param(params, IGBPs), gof, output=df_out) +end + + +## 一个站点的率定 +""" + ModelCalib(df::AbstractDataFrame, par0::AbstractETParam; + IGBPcode=nothing, maxn=2500, of_gof=:NSE, kw...) +""" +function ModelCalib(df::AbstractDataFrame, par0::AbstractETParam, parNames; + IGBPcode=nothing, maxn=2500, of_gof=:NSE, kw...) + lower, upper = get_bounds(parNames) + x0 = select_param(par0, parNames) + + sceua(theta -> -model_goal(df, theta, parNames; IGBPcode, of_gof, kw...), + x0, lower, upper; maxn, kw...) # theta, goal, flag +end + + +## gof of ET and GPP for one IGBP +# cal df_out by Forcing and Param +function model_goal(df, theta, parNames; IGBPcode=nothing, of_gof=:NSE, verbose=false) + par = theta2par(theta, parNames) + :hc ∉ parNames && !isnothing(IGBPcode) && (par.hc = hc_raw[IGBPcode]) + + dobs = df[!, [:GPP_obs, :ET_obs]] + dsim = PMLV2_sites(df; par) + # dsim = PMLV2(df; par) # the exact MATLAB version + + ## 8-day, yearly, yearly_anom + info_GPP = GOF(dobs.GPP_obs, dsim.GPP) + info_ET = GOF(dobs.ET_obs, dsim.ET) + + if verbose + indexes = [:KGE, :NSE, :R2, :RMSE, :bias, :bias_perc] + info_GPP = info_GPP[indexes] |> round2 + info_ET = info_ET[indexes] |> round2 + end + gof = [info_ET[of_gof], info_GPP[of_gof]] + weighted_mean(gof, [1.0, 1.0]) # goal +end + + +## gof of ET and GPP for ALL IGBP +function model_gof(df_out::DataFrame; all=true) + fun_et(d) = GOF(d.ET_obs, d.ET) + fun_gpp(d) = GOF(d.GPP_obs, d.GPP) + + lst = groupby(df_out, :IGBP) + + res = map_df_tuple(fun_et, lst) + all && push!(res, (; IGBP="ALL", fun_et(df_out)...)) + gof_et = DataFrame(res) + + res = map_df_tuple(fun_gpp, lst) + all && push!(res, (; IGBP="ALL", fun_gpp(df_out)...)) + gof_gpp = DataFrame(res) + (; ET=gof_et, GPP=gof_gpp) +end diff --git a/src/Optim/ModelCalib.jl b/src/Optim/ModelCalib.jl deleted file mode 100644 index b418d6e..0000000 --- a/src/Optim/ModelCalib.jl +++ /dev/null @@ -1,54 +0,0 @@ -export ModelCalib, ModelCalib_IGBPs -export fread, fwrite, melt_list - -using RTableTools -import Ipaper: par_map -include("DataType.jl") -include("model_gof.jl") -# include("PMLV2_sites.jl") - - -function ModelCalib_IGBPs(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"] - IGBPs = unique(data.IGBPname) |> sort - - @time params = par_map(IGBP -> begin - df = data[data.IGBP.==IGBP, vars] - IGBPcode = df.IGBPcode[1] - theta, _, _ = ModelCalib(df, par0; IGBPcode, of_gof, maxn, verbose=false) - theta - end, IGBPs; progress=false) - - # printstyled("[i=$i, IGBP = $IGBP] \n", bold=true, color=:green, underline=false) - res = map(i -> begin - IGBP = IGBPs[i] - theta = params[i] - df = data[data.IGBP.==IGBP, vars] - model_goal(df, theta; verbose=true) - - par = theta2par(theta) - r = PMLV2_sites(df; par) - cbind(df[:, [:site, :date, :ET_obs, :GPP_obs]], r) - end, eachindex(params)) - df_out = melt_list(res; IGBP=IGBPs) - - gof = model_gof(df_out) - (; param=theta2param(params, IGBPs), gof, output=df_out) -end - - -## 一个站点的率定 -""" - ModelCalib(df::AbstractDataFrame, par0::AbstractETParam; - IGBPcode=nothing, maxn=2500, of_gof=:NSE, kw...) -""" -function ModelCalib(df::AbstractDataFrame, par0::AbstractETParam; IGBPcode=nothing, maxn=2500, of_gof=:NSE, kw...) - parRanges = get_bounds(par0) - lower = parRanges[:, 1] - upper = parRanges[:, 2] - theta0 = collect(par0) # must be vector - - sceua(theta -> -model_goal(df, theta; IGBPcode, of_gof, kw...), - theta0, lower, upper; maxn, kw...) # theta, goal, flag -end diff --git a/src/Optim/model_gof.jl b/src/Optim/model_gof.jl deleted file mode 100644 index 2d1683f..0000000 --- a/src/Optim/model_gof.jl +++ /dev/null @@ -1,56 +0,0 @@ -function map_df_tuple(fun::Function, lst::GroupedDataFrame{DataFrame}, args...; kw...) - n = length(lst) - _keys = keys(lst) - map(i -> begin - d = lst[i] - key = NamedTuple(_keys[i]) - r = fun(d, args...; kw...) - (; key..., r...) - end, 1:n) -end - - -## 统计每种植被类型的模拟效果 -function model_gof(df_out::DataFrame; all=true) - fun_et(d) = GOF(d.ET_obs, d.ET) - fun_gpp(d) = GOF(d.GPP_obs, d.GPP) - - lst = groupby(df_out, :IGBP) - - res = map_df_tuple(fun_et, lst) - all && push!(res, (; IGBP="ALL", fun_et(df_out)...)) - gof_et = DataFrame(res) - - res = map_df_tuple(fun_gpp, lst) - all && push!(res, (; IGBP="ALL", fun_gpp(df_out)...)) - gof_gpp = DataFrame(res) - (; ET=gof_et, GPP=gof_gpp) -end - - -function model_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 = Param_PMLV2(theta...) - - dobs = df[!, [:GPP_obs, :ET_obs]] - # dsim = PMLV2(df; par) # the exact MATLAB version - dsim = PMLV2_sites(df; par) - - ## 8-day, yearly, yearly_anom - info_GPP = GOF(dobs.GPP_obs, dsim.GPP) - info_ET = GOF(dobs.ET_obs, dsim.ET) - - if verbose - indexes = [:KGE, :NSE, :R2, :RMSE, :bias, :bias_perc] - info_GPP = info_GPP[indexes] |> round2 - info_ET = info_ET[indexes] |> round2 - end - gof = [info_ET[of_gof], info_GPP[of_gof]] - goal = weighted_mean(gof, [1, 0.5]) - goal -end - - -export model_goal -export model_gof, map_df_tuple diff --git a/src/PML.jl b/src/PML.jl index e6c942b..e6c7a2b 100644 --- a/src/PML.jl +++ b/src/PML.jl @@ -3,16 +3,22 @@ module PML export PMLV2, PMLV2_sites, photosynthesis, cal_Ei_Dijk2021, T_adjust_Vm25, f_VPD_Zhang2019 -export DataFrame, GOF + export file_FLUXNET_CRO, file_FLUXNET_CRO_USTwt +export DataFrame, GOF +export fread, fwrite, melt_list + using Parameters using FieldMetadata using DataFrames using Statistics -import HydroTools: cal_Uz, Cp, atm, GOF, sceua +using RTableTools using DocStringExtensions +import HydroTools: cal_Uz, Cp, atm, GOF, sceua +import Ipaper: par_map + ## global data dir_proj = "$(@__DIR__)/.." file_FLUXNET_CRO = "$dir_proj/data/CRO/FLUXNET_CRO" |> abspath @@ -21,18 +27,16 @@ file_FLUXNET_CRO_USTwt = "$dir_proj/data/CRO/FLUXNET_CRO_US-Twt" |> abspath include("main_Ipaper.jl") include("Params.jl") -include("Optim/ModelCalib.jl") +include("ModelCalib.jl") include("ET_helper.jl") +include("water_constrain.jl") +include("photosynthesis.jl") # include("PET_equilibrium.jl") # include("Ei_EvapIntercepted.jl") # include("Ec_CanopyTrans.jl") # include("Es_EvapSoil.jl") -include("water_constrain.jl") -include("photosynthesis.jl") # include("model_PMLV2.jl") - - """ PMLV2 (Penman–Monteith–Leuning Version 2) Evapotranspiration model diff --git a/src/Params.jl b/src/Params.jl index f532896..f8fc31f 100644 --- a/src/Params.jl +++ b/src/Params.jl @@ -58,26 +58,30 @@ function Base.collect(par::AbstractETParam) [getfield(par, f) for f in fieldnames(typeof(par))] end -function get_bounds(par::AbstractETParam, parNames::Vector{Symbol}=ParNames) +function get_bounds(parNames::Vector{Symbol}=ParNames) inds = match2(parNames, ParNames).I_y # 选择的参数 - bs = bounds(par)[inds] - x0 = vcat([getfield(par, f) for f in parNames]...) + bs = bounds(Param_PMLV2)[inds] lower = vcat(map(x -> x[1], bs)...) upper = vcat(map(x -> x[2], bs)...) - x0, lower, upper + lower, upper end # theta2par(theta) = Param_PMLV2(theta...) -function theta2par(theta::Vector, par::AbstractETParam, parNames::Vector{Symbol}) +function theta2par!(theta::Vector, par::AbstractETParam, parNames::Vector{Symbol}) k = 0 for key in parNames x = getfield(par, key) n = length(x) - inds = k+1:k+n + inds = n == 1 ? k+1 : k+1:k+n k += n setfield!(par, key, theta[inds]) end + return par +end + +function theta2par(theta::Vector, parNames::Vector{Symbol}) + theta2par!(theta, deepcopy(par0), parNames) end function select_param(par::AbstractETParam, parNames) @@ -88,11 +92,8 @@ end hc_raw = [10, 10, 10, 10, 10, 1, 1, 5, 5, 0.2, 1, 0.5, 10, 1, 0.01, 0.05, 0.01] FT = Float64 - par0 = Param_PMLV2() theta0 = collect(par0) -# parNames = fieldnames(Param_PMLV2) -# parRanges = get_bounds(par0) function theta2param(params::Vector{Vector{T}}, IGBPs) where {T<:Real} parNames = fieldnames(Param_PMLV2) |> collect @@ -103,7 +104,7 @@ function theta2param(params::Vector{Vector{T}}, IGBPs) where {T<:Real} end -export bounds, Param_PMLV2 -export parRanges, parNames, theta0, par0, hc_raw -export get_bounds, select_param +export bounds, get_bounds, select_param +export Param_PMLV2, ParNames +export theta0, par0, hc_raw export param_PML, theta2par, theta2param diff --git a/test/test-PMLV2.jl b/test/test-PMLV2.jl index abb862f..2788507 100644 --- a/test/test-PMLV2.jl +++ b/test/test-PMLV2.jl @@ -1,8 +1,9 @@ using PML, Test, Ipaper, DataFrames df_out, df, _par = deserialize(file_FLUXNET_CRO_USTwt) -par = Param_PMLV2(;_par..., hc=0.5) - +df.GPP_obs = df.GPPobs +df.ET_obs = df.ETobs +par = Param_PMLV2(; _par..., hc=0.5) r = PMLV2_sites(df; par) @testset "PMLV2 scalar" begin @@ -32,18 +33,21 @@ end end -# @testset "ModelCalib" begin -# df.GPP_obs = df.GPPobs -# df.ET_obs = df.ETobs -# @time _theta, goal, flag = ModelCalib(df, par0; IGBPcode=df.IGBPcode[1], maxn=2500) -# goal = model_goal(df, _theta; verbose=true) -# @test goal > 0.55 # mean(KGE_GPP, KGE_ET) -# end - -# @testset "model_gof" begin -# d = DataFrame(; IGBP=["CRO", "CRO", "CRO"], -# ET_obs=[1.0, 2, 3], ET=[0.5, 0.6, 0.7], -# GPP_obs=[1.0, 2, 3], GPP=[1, 1.5, 2.0]) -# @test size(model_gof(d).ET, 1) == 2 -# @test size(model_gof(d; all=false).ET, 1) == 1 -# end +@testset "model_gof" begin + d = DataFrame(; IGBP=["CRO", "CRO", "CRO"], + ET_obs=[1.0, 2, 3], ET=[0.5, 0.6, 0.7], + GPP_obs=[1.0, 2, 3], GPP=[1, 1.5, 2.0]) + @test size(model_gof(d).ET, 1) == 2 + @test size(model_gof(d; all=false).ET, 1) == 1 +end + + +@testset "ModelCalib" begin + parNames = [ + :α, :η, :g1, :Am_25, :VPDmin, :VPDmax, :D0, :kQ, :kA, :S_sls, :fER0 # :hc + ] + + @time _theta, goal, flag = ModelCalib(df, par0, parNames; IGBPcode=df.IGBPcode[1], maxn=2500) + goal = model_goal(df, _theta, parNames; verbose=true) + @test goal > 0.55 # mean(KGE_GPP, KGE_ET) +end