From 9db88ca99434bedd0ffd7d974f313df201cbc606 Mon Sep 17 00:00:00 2001 From: lalonso Date: Fri, 11 Jul 2025 12:59:36 +0100 Subject: [PATCH 01/27] fixes lts support and add CI check for it --- .github/workflows/CI.yml | 9 +++++++++ Project.toml | 4 ++-- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 41c35edbf..202344a1e 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -15,11 +15,20 @@ jobs: fail-fast: false matrix: version: + - '1.10' - '1.11' os: - ubuntu-latest + - macOS-latest + - windows-latest arch: - x64 + allow_failure: [false] + include: + - version: 'nightly' + os: ubuntu-latest + arch: x64 + allow_failure: true steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 diff --git a/Project.toml b/Project.toml index 99c0c43f7..d94d68016 100644 --- a/Project.toml +++ b/Project.toml @@ -25,12 +25,12 @@ TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9" [compat] Accessors = "0.1.42" Crayons = "4.1.1" -Dates = "1.11.0" +Dates = "1.10, 1.11.0" NaNStatistics = "0.6.50" Pkg = "1.10.0" Reexport = "1.2.2" TypedTables = "1.4.6" -julia = "1.9" +julia = "1.10" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 4d6dc3252fa2a079bac61bd8d0b4bbccc0a6a882 Mon Sep 17 00:00:00 2001 From: lalonso Date: Fri, 11 Jul 2025 13:00:44 +0100 Subject: [PATCH 02/27] fixes path name --- lib/SindbadSetup/src/setupExperimentInfo.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/lib/SindbadSetup/src/setupExperimentInfo.jl b/lib/SindbadSetup/src/setupExperimentInfo.jl index 6015da3b3..340e13ede 100644 --- a/lib/SindbadSetup/src/setupExperimentInfo.jl +++ b/lib/SindbadSetup/src/setupExperimentInfo.jl @@ -104,8 +104,7 @@ Saves or skips saving the experiment configuration to a file. function saveInfo end function saveInfo(info, ::DoSaveInfo) - info_path = joinpath(info.output.dirs.settings, "info.jld - 2") + info_path = joinpath(info.output.dirs.settings, "info.jld2") showInfo(saveInfo, @__FILE__, @__LINE__, "saving info to `$(info_path)`") @save info_path info return nothing From 6972dc6ab5cc4919730ca1b20ed0ddb02c346a4d Mon Sep 17 00:00:00 2001 From: lalonso Date: Fri, 11 Jul 2025 13:06:41 +0100 Subject: [PATCH 03/27] in windows, create path if it does not exists --- lib/SindbadSetup/src/setupOutput.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/SindbadSetup/src/setupOutput.jl b/lib/SindbadSetup/src/setupOutput.jl index b0553986f..c9c9cf034 100644 --- a/lib/SindbadSetup/src/setupOutput.jl +++ b/lib/SindbadSetup/src/setupOutput.jl @@ -255,6 +255,7 @@ Saves a copy of the experiment settings to the output folder. function saveExperimentSettings(info) sindbad_experiment = info.temp.experiment.dirs.sindbad_experiment showInfo(saveExperimentSettings, @__FILE__, @__LINE__, "saving Experiment JSON Settings to : $(info.output.dirs.settings)") + mkpath(dirname(joinpath(info.output.dirs.settings, split(sindbad_experiment, path_separator)[end]))) # create output directory if it does not exist cp(sindbad_experiment, joinpath(info.output.dirs.settings, split(sindbad_experiment, path_separator)[end]); force=true) From 23956ac080514e9b8f0a5eaf274ff9f6d5ac51e3 Mon Sep 17 00:00:00 2001 From: lalonso Date: Fri, 11 Jul 2025 13:12:01 +0100 Subject: [PATCH 04/27] don't include Sindbad.Models in full_name --- lib/SindbadSetup/src/setupParameters.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/SindbadSetup/src/setupParameters.jl b/lib/SindbadSetup/src/setupParameters.jl index b4edfc7f3..d817a4f66 100644 --- a/lib/SindbadSetup/src/setupParameters.jl +++ b/lib/SindbadSetup/src/setupParameters.jl @@ -75,7 +75,8 @@ function getParameters(selected_models::Tuple, num_type, model_timestep; return_ upper = [constrains[i][2] for i in 1:nbounds] model = [Symbol(supertype(getproperty(Models, m))) for m in model_approach] - name_full = [join((model[i], name[i]), ".") for i in 1:nbounds] + model_str = string.(model) + name_full = [join((last(split(model_str[i], ".")), name[i]), ".") for i in 1:nbounds] approach_func = [getfield(Models, m) for m in model_approach] model_prev = model_approach[1] m_id = findall(x-> x==model_prev, model_names_list)[1] From c978c18c0c842971fda684fc36edf3d59a1a53a3 Mon Sep 17 00:00:00 2001 From: lalonso Date: Fri, 11 Jul 2025 13:33:10 +0100 Subject: [PATCH 05/27] info.temp missing field fix --- lib/SindbadSetup/src/setupHybridML.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SindbadSetup/src/setupHybridML.jl b/lib/SindbadSetup/src/setupHybridML.jl index d9c1b02c2..13f85c6f4 100644 --- a/lib/SindbadSetup/src/setupHybridML.jl +++ b/lib/SindbadSetup/src/setupHybridML.jl @@ -152,7 +152,7 @@ function setHybridInfo(info::NamedTuple) info = setTupleSubfield(info, :hybrid, (:replace_value_for_gradient, info.temp.helpers.numbers.num_type(replace_value_for_gradient))) - covariates_path = getAbsDataPath(info, info.settings.hybrid.covariates.path) + covariates_path = getAbsDataPath(info.temp, info.settings.hybrid.covariates.path) covariates = (; path=covariates_path, variables=info.settings.hybrid.covariates.variables) info = setTupleSubfield(info, :hybrid, (:covariates, covariates)) info = setTupleSubfield(info, :hybrid, (:random_seed, info.settings.hybrid.random_seed)) From e431b2c61cfe7682e6d2d812cb9fb8fc771315f1 Mon Sep 17 00:00:00 2001 From: lalonso Date: Tue, 15 Jul 2025 13:32:54 +0100 Subject: [PATCH 06/27] fix destination path for experiment info --- lib/SindbadSetup/src/setupOutput.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/SindbadSetup/src/setupOutput.jl b/lib/SindbadSetup/src/setupOutput.jl index c9c9cf034..29aee30d3 100644 --- a/lib/SindbadSetup/src/setupOutput.jl +++ b/lib/SindbadSetup/src/setupOutput.jl @@ -255,9 +255,9 @@ Saves a copy of the experiment settings to the output folder. function saveExperimentSettings(info) sindbad_experiment = info.temp.experiment.dirs.sindbad_experiment showInfo(saveExperimentSettings, @__FILE__, @__LINE__, "saving Experiment JSON Settings to : $(info.output.dirs.settings)") - mkpath(dirname(joinpath(info.output.dirs.settings, split(sindbad_experiment, path_separator)[end]))) # create output directory if it does not exist + destination_experiment_file = last(split(last(split(sindbad_experiment, path_separator)), "/")) cp(sindbad_experiment, - joinpath(info.output.dirs.settings, split(sindbad_experiment, path_separator)[end]); + joinpath(info.output.dirs.settings, destination_experiment_file); force=true) for k ∈ keys(info.settings.experiment.basics.config_files) v = getfield(info.settings.experiment.basics.config_files, k) From c7817eec2d8d8335e47a3e7f0b5430434f3551d9 Mon Sep 17 00:00:00 2001 From: lazarusA Date: Fri, 18 Jul 2025 10:29:18 +0200 Subject: [PATCH 07/27] get first global parameters --- lib/SindbadML/src/mapParams.jl | 137 +++++++++++++++++++++++++++++++++ 1 file changed, 137 insertions(+) create mode 100644 lib/SindbadML/src/mapParams.jl diff --git a/lib/SindbadML/src/mapParams.jl b/lib/SindbadML/src/mapParams.jl new file mode 100644 index 000000000..a942e74d2 --- /dev/null +++ b/lib/SindbadML/src/mapParams.jl @@ -0,0 +1,137 @@ +export stackFeatures +export mapParamsPFT +export mapParamsAll + +""" + mapNNYaxPFT(out_ps, in_pft; trainedNN, lower_bound, upper_bound) + +`mapCube` wrapper function that outpus `out_ps` new parameters. + +- `in_pft`: input Plant Functional Type (PFT) cube entry. +- `trainedNN`: A trained neural network. +- `lower_bound`: parameters' lower bounds. +- `upper_bound`: parameters' upper bounds. +""" +function mapNNYaxPFT(out_ps, in_pft; trainedNN, lower_bound, upper_bound) + x_in = Flux.onehot(only(in_pft), 1:17, 17) + new_ps = scaleToBounds.(trainedNN(x_in), lower_bound, upper_bound) + out_ps[:] = new_ps # ? doing explicitly `.=` also works! +end + +""" + mapNNYaxAll(out_ps, in_pft, in_kg, in_all_args; trainedNN, lower_bound, upper_bound) + +`mapCube` wrapper function that outputs `out_ps` new parameters using all covariates. + +Arguments: +- `out_ps`: output parameters. +- `in_pft`: input Plant Functional Type (PFT) cube entry. +- `in_kg`: input Koeppen-Geiger (KG) cube entry. +- `in_all_args`: all other input arguments. +- `trainedNN`: A trained neural network. +- `lower_bound`: parameters' lower bounds. +- `upper_bound`: parameters' upper bounds. +""" +function mapNNYaxAll(out_ps, in_pft, in_kg, in_all_args; trainedNN, lower_bound, upper_bound) + x_in_pft = Flux.onehot(only(in_pft), 1:17, 17) + x_in_kg = Flux.onehot(only(in_kg), 1:32, 32) + x_in = reduce(vcat, [x_in_kg, x_in_pft, in_all_args...]) + # @show length(x_in) + new_ps = scaleToBounds.(trainedNN(x_in), lower_bound, upper_bound) + out_ps[:] = new_ps # ? doing explicitly `.=` also works! +end + +""" + stackFeatures(pft, kg, add_args...; up_bound=17, veg_cat=false, clim_cat=false) + +Stack all features into a vector. + +- `pft`: LandCover type (Plant Functional Type). Any entry not in 1:17 would be set to the last index, this includes NaN! Last index is water/NaN +- `kg`: Koeppen-Geiger climater type. Any entry not in 1:32 would be set to the last index, this includes NaN! Last index is water/NaN +- `add_args`: all non-categorical additional features +- `up_bound`: last index class, the range goes from `1:up_bound`, and any case not in that range uses the `up_bound` value. For `PFT` use `17`. +- `veg_cat`: `true` or `false`. + +Returns a vector. +""" +function stackFeatures(pft, kg, add_args...; up_bound=17, veg_cat=false, clim_cat=false) + veg_onehot = oneHotPFT(pft, up_bound, veg_cat) + if !clim_cat + return reduce(vcat, [veg_onehot, add_args...]) + else + kg_onehot = Flux.onehot(kg, 1:32, 32) + return reduce(vcat, [kg_onehot, veg_onehot, add_args...]) + end +end + + +""" + mapParamsAll(incubes, trainedNN, lower_bound, upper_bound, ps_names, path; metadata_global = Dict()) + +Compute all parameters using a neural network and all input covariates. + +Arguments: +- `incubes`: input covariates cubes. Firsts ones should be `pft` and `kc`. +- `trainedNN`: A trained neural network. +- `lower_bound`: parameters' lower bounds. +- `upper_bound`: parameters' upper bounds. +- `ps_names`: parameter names +- `path`: output path for the cube. +- `metadata_global`: global metadata to be added to the output cube. +""" +function mapParamsAll(incubes, trainedNN, lower_bound, upper_bound, ps_names, path; metadata_global = Dict()) + + indims = [InDims(; filter=AllNaN()), InDims(; filter=AllNaN()), InDims("Variables"; filter=AllNaN())] + + properties = Dict{String, Any}( + "name"=> "parameters", + "description" => "neural network spatial parameters' estimations" + ) + properties = merge(properties, metadata_global) + mapCube(mapNNYaxAll, (incubes...,); # ! additional input function arguments + trainedNN, + lower_bound, + upper_bound, + indims = indims, + outdims = OutDims(Dim{:parameter}(String.(ps_names)); + path=path, + outtype=Float32, + properties, + overwrite=true), + ) +end + +""" + mapParamsPFT(inPFTcube, trainedNN, lower_bound, upper_bound, ps_names, path; metadata_global = Dict()) + +Compute parameters using a neural network and the PFT covariates. + +Arguments: +- `inPFTcube`: input Plant Functional Type (PFT) cube. +- `trainedNN`: A trained neural network. +- `lower_bound`: parameters' lower bounds. +- `upper_bound`: parameters' upper bounds. +- `ps_names`: parameter names. +- `path`: output path for the cube. +- `metadata_global`: global metadata to be added to the output cube. +""" +function mapParamsPFT(inPFTcube, trainedNN, lower_bound, upper_bound, ps_names, path; metadata_global = Dict()) + indims = [InDims(; filter=AllNaN())] + + properties = Dict{String, Any}( + "name"=> "parameters", + "description" => "neural network spatial parameters' estimations" + ) + properties = merge(properties, metadata_global) + mapCube(mapNNYaxPFT, (inPFTcube,); # ! additional input function arguments + trainedNN, + lower_bound, + upper_bound, + indims = indims, + outdims = OutDims(Dim{:parameter}(String.(ps_names)); + path=path, + outtype=Float32, + properties, + overwrite=true), + ) +end \ No newline at end of file From 08ae567883937f4d2169b6d882c537c2357e4158 Mon Sep 17 00:00:00 2001 From: lazarusA Date: Fri, 18 Jul 2025 10:36:38 +0200 Subject: [PATCH 08/27] and now do independent model runs for each pixel --- lib/SindbadTEM/src/runTEMCube.jl | 102 +++++++++++++++++++++++++++---- 1 file changed, 90 insertions(+), 12 deletions(-) diff --git a/lib/SindbadTEM/src/runTEMCube.jl b/lib/SindbadTEM/src/runTEMCube.jl index c46a57fb4..db72fd4e5 100644 --- a/lib/SindbadTEM/src/runTEMCube.jl +++ b/lib/SindbadTEM/src/runTEMCube.jl @@ -17,13 +17,11 @@ run the SINBAD CORETEM for a given location function coreTEMYax(selected_models, loc_forcing, loc_land, tem_info) loc_forcing_t = getForcingForTimeStep(loc_forcing, deepcopy(loc_forcing), 1, tem_info.vals.forcing_types) - - spinup_forcing = getAllSpinupForcing(loc_forcing, tem_info.spinup_sequence, tem_info.model_helpers); + spinup_forcing = getAllSpinupForcing(loc_forcing, tem_info.spinup_sequence, tem_info); land_prec = definePrecomputeTEM(selected_models, loc_forcing_t, loc_land, tem_info.model_helpers) - - land_spin = spinupTEM(selected_models, spinup_forcing, loc_forcing_t, land_prec, tem_info, run_helpers.tem_info.run.spinup_TEM) - + land_prec = precomputeTEM(selected_models, loc_forcing_t, land_prec, tem_info.model_helpers) + land_spin = spinupTEM(selected_models, spinup_forcing, loc_forcing_t, land_prec, tem_info) # ? not good? see git diff # land_spin = land_prec land_time_series = timeLoopTEM(selected_models, loc_forcing, loc_forcing_t, land_spin, tem_info, tem_info.run.debug_model) @@ -43,14 +41,22 @@ end - `selected_models`: a tuple of all models selected in the given model structure - `forcing_vars`: forcing variables """ -function TEMYax(map_cubes...;selected_models::Tuple, forcing_vars::AbstractArray, loc_land::NamedTuple, output_vars, tem::NamedTuple) +function TEMYax(map_cubes...;selected_models::Tuple, forcing_vars, loc_land::NamedTuple, output_vars, tem::NamedTuple, clean_data) + outputs, inputs = unpackYaxForward(map_cubes; output_vars, forcing_vars) + + # ? apply clean_data fields to input data points + _data_fill, _forcing_default_info, _num_type, _forcing_vars_info = clean_data + + inputs = map(enumerate(inputs)) do (i, in_forcing) + _data_info = getCombinedNamedTuple(_forcing_default_info, _forcing_vars_info[i]) + map(data_point -> cleanData(data_point, _data_fill, _data_info, _num_type), in_forcing) + end loc_forcing = (; Pair.(forcing_vars, inputs)...) land_out = coreTEMYax(selected_models, loc_forcing, loc_land, tem) - i = 1 foreach(output_vars) do var_pair - # @show i, var_pair + # @show i, var_pair, size(outputs[i]) data = land_out[first(var_pair)][last(var_pair)] fillOutputYax(outputs[i], data) i += 1 @@ -79,20 +85,92 @@ function runTEMYax(selected_models::Tuple, forcing::NamedTuple, info::NamedTuple run_helpers = prepTEM(forcing, info); loc_land = deepcopy(run_helpers.loc_land); + _data_fill = 0.0f0 + _forcing_default_info = info.settings.forcing.default_forcing + _num_type = Val{info.helpers.numbers.num_type}() + _forcing_vars_info = info.settings.forcing.variables + outcubes = mapCube(TEMYax, (incubes...,); selected_models=selected_models, forcing_vars=forcing.variables, - output_vars = run_helpers.output_vars, + output_vars=run_helpers.output_vars, loc_land=loc_land, tem=run_helpers.tem_info, + clean_data=(; _data_fill, _forcing_default_info, _num_type, _forcing_vars_info), indims=indims, outdims=run_helpers.output_dims, max_cache=info.settings.experiment.exe_rules.yax_max_cache, - ispar=true) + ispar=false, + ) return outcubes end +""" + psTEMYax(in_pixel_cube...; selected_models::Tuple, param_to_index, forcing_vars, loc_land::NamedTuple, output_vars, tem::NamedTuple, clean_data) +""" +function psTEMYax(in_pixel_cube...; selected_models::Tuple, param_to_index, forcing_vars, loc_land::NamedTuple, + output_vars, tem::NamedTuple, clean_data) + + outputs, inputs = unpackYaxForward(in_pixel_cube[1:end-1]; output_vars, forcing_vars) + in_new_params = last(in_pixel_cube) + # ? apply clean_data fields to input data points + _data_fill, _forcing_default_info, _num_type, _forcing_vars_info = clean_data + + inputs = map(enumerate(inputs)) do (i, in_forcing) + _data_info = getCombinedNamedTuple(_forcing_default_info, _forcing_vars_info[i]) + map(data_point -> cleanData(data_point, _data_fill, _data_info, _num_type), in_forcing) + end + loc_forcing = (; Pair.(forcing_vars, inputs)...) + updated_models = updateModelParameters(param_to_index, selected_models, in_new_params) + land_out = coreTEMYax(updated_models, loc_forcing, loc_land, tem) + i = 1 + foreach(output_vars) do var_pair + # @show i, var_pair, size(outputs[i]) + data = land_out[first(var_pair)][last(var_pair)] + fillOutputYax(outputs[i], data) + i += 1 + end +end + + +""" + run_psTEMYax(selected_models::Tuple, forcing::NamedTuple, in_cube_params, tbl_params, info::NamedTuple) + +""" +function run_psTEMYax(selected_models::Tuple, forcing::NamedTuple, in_cube_params, tbl_params, info::NamedTuple) + + # forcing/input information + in_cubes_forcing = forcing.data + in_cubes_all = (in_cubes_forcing..., in_cube_params) + indims = forcing.dims + indims = (indims..., InDims((YAXArrays.ByName("parameter"),), Array, (AllNaN(),), missing)) + param_to_index = getParameterIndices(selected_models, tbl_params) + # information for running model + run_helpers = prepTEM(forcing, info) + loc_land = deepcopy(run_helpers.loc_land) + + _data_fill = 0.0f0 + _forcing_default_info = info.settings.forcing.default_forcing + _num_type = Val{info.helpers.numbers.num_type}() + _forcing_vars_info = info.settings.forcing.variables + + outcubes = mapCube(psTEMYax, + (in_cubes_all...,); + selected_models=selected_models, + param_to_index=param_to_index, + forcing_vars=forcing.variables, + output_vars=run_helpers.output_vars, + loc_land=loc_land, + tem=run_helpers.tem_info, + clean_data=(; _data_fill, _forcing_default_info, _num_type, _forcing_vars_info), + indims=indims, + outdims=run_helpers.output_dims, + max_cache=info.settings.experiment.exe_rules.yax_max_cache, + ispar=true, + ) + return outcubes +end """ unpackYaxForward(args; tem::NamedTuple, forcing_vars::AbstractArray) @@ -124,13 +202,13 @@ fills the output array position with the input data/vector - `xin`: input data/vector """ function fillOutputYax(xout, xin) - if ndims(xout) == ndims(xin) + if ndims(xout) == ndims(xin) && length(xin[1]) == 1 for i ∈ eachindex(xin) xout[i] = xin[i][1] end else for i ∈ CartesianIndices(xin) - xout[:, i] .= xin[i] + xout[:, i] = xin[i] end end end \ No newline at end of file From d98225863d25e211c4f26eb44f04bbef7a42c640 Mon Sep 17 00:00:00 2001 From: lazarusA Date: Fri, 18 Jul 2025 10:50:59 +0200 Subject: [PATCH 09/27] only clean data if the output type is an array, otherwise, the cleaning its done at loading, otherwise we will run out of memory for global datasets --- lib/SindbadData/src/utilsData.jl | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/lib/SindbadData/src/utilsData.jl b/lib/SindbadData/src/utilsData.jl index eb4cdcafc..5ea778a8c 100644 --- a/lib/SindbadData/src/utilsData.jl +++ b/lib/SindbadData/src/utilsData.jl @@ -430,16 +430,18 @@ function subsetAndProcessYax(yax, forcing_mask, tar_dims, _data_info, info, ::Va end #todo mean of the data instead of zero or nan - vfill = 0.0 - if fill_nan - vfill = NaN - end - vNT = Val{num_type}() - if clean_data - yax = mapCleanData(yax, yax_qc, vfill, bounds_qc, _data_info, vNT) - else - yax = map(yax_point -> replaceInvalid(yax_point, vfill), yax) - yax = num_type.(yax) + if info.settings.experiment.model_output.output_array_type == "array" + vfill = 0.0 + if fill_nan + vfill = NaN + end + vNT = Val{num_type}() + if clean_data + yax = mapCleanData(yax, yax_qc, vfill, bounds_qc, _data_info, vNT) + else + yax = map(yax_point -> replaceInvalid(yax_point, vfill), yax) + yax = num_type.(yax) + end end return yax end From 0cf26f4a6230279963591a0a447d785168a12d53 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Mon, 21 Jul 2025 11:53:35 +0200 Subject: [PATCH 10/27] test global run --- examples/exp_global_runs/exp_global_runs.jl | 1 + .../settings_global_runs/experiment.json | 89 ++++++ .../settings_global_runs/forcing.json | 184 +++++++++++++ .../forcing_f_backup.json | 184 +++++++++++++ .../settings_global_runs/model_structure.json | 258 ++++++++++++++++++ 5 files changed, 716 insertions(+) create mode 100644 examples/exp_global_runs/exp_global_runs.jl create mode 100644 examples/exp_global_runs/settings_global_runs/experiment.json create mode 100644 examples/exp_global_runs/settings_global_runs/forcing.json create mode 100644 examples/exp_global_runs/settings_global_runs/forcing_f_backup.json create mode 100644 examples/exp_global_runs/settings_global_runs/model_structure.json diff --git a/examples/exp_global_runs/exp_global_runs.jl b/examples/exp_global_runs/exp_global_runs.jl new file mode 100644 index 000000000..0c1faa85b --- /dev/null +++ b/examples/exp_global_runs/exp_global_runs.jl @@ -0,0 +1 @@ +2 + 2 \ No newline at end of file diff --git a/examples/exp_global_runs/settings_global_runs/experiment.json b/examples/exp_global_runs/settings_global_runs/experiment.json new file mode 100644 index 000000000..07e9880b4 --- /dev/null +++ b/examples/exp_global_runs/settings_global_runs/experiment.json @@ -0,0 +1,89 @@ +{ + "basics": { + "config_files": { + "forcing": "forcing.json", + "model_structure": "model_structure.json" + }, + "domain": "Global", + "name": "Fix_ALL_ALL", + "time": { + "date_begin": "2002-01-01", + "date_end": "2022-12-30", + "temporal_resolution": "day", + "timesteps_in_day": 1, + "timesteps_in_year": 365 + } + }, + "exe_rules": { + "input_array_type": "yax_array", + "input_data_backend": "zarr", + "land_output_type": "YAXArray", + "longtuple_size": 10, + "model_array_type": "static_array", + "model_number_type": "Float32", + "parallelization": "threads", + "tolerance": 1.0e-2, + "yax_max_cache": 2e9 + }, + "flags": { + "calc_cost": false, + "catch_model_errors": false, + "debug_model": false, + "filter_nan_pixels": false, + "inline_update": false, + "run_forward": true, + "run_optimization": false, + "save_info": true, + "spinup_TEM": true, + "store_spinup": false, + "use_forward_diff": false + }, + "model_output": { + "depth_dimensions": { + "d_cEco": 8, + "d_snow": "snowW", + "d_soil": "soilW", + "d_tws": "TWS" + }, + "format": "zarr", + "output_array_type": "YAXArray", + "path": null, + "save_single_file": true, + "variables": { + "diagnostics.c_eco_k_f_soilT": null, + "diagnostics.c_eco_k_f_soilW": "d_cEco", + "fluxes.auto_respiration": null, + "fluxes.evapotranspiration": null, + "fluxes.gpp": null, + "fluxes.nee": null, + "fluxes.transpiration": null, + "pools.cEco": "d_cEco", + "diagnostics.c_Fire_k": "d_cEco", + "diagnostics.c_Fire_cci": "d_cEco", + "diagnostics.c_Veg_Mortality": "d_cEco", + "diagnostics.c_Fire_Flux": "d_cEco", + "states.fAPAR_bare": null + } + }, + "model_spinup": { + "restart_file": null, + "sequence": [ + { + "forcing": "all_years", + "n_repeat": 1, + "spinup_mode": "all_forward_models" + }, + { + "forcing": "day_MSC", + "n_repeat": 1, + "spinup_mode": "all_forward_models" + }, + { + "forcing": "first_year", + "n_repeat": 200, + "spinup_mode": "all_forward_models" + } + ] + } + } + \ No newline at end of file diff --git a/examples/exp_global_runs/settings_global_runs/forcing.json b/examples/exp_global_runs/settings_global_runs/forcing.json new file mode 100644 index 000000000..5f174c075 --- /dev/null +++ b/examples/exp_global_runs/settings_global_runs/forcing.json @@ -0,0 +1,184 @@ +{ + "data_dimension": { + "time": "Ti", + "permute": null, + "space": ["longitude", "latitude"] + }, + "default_forcing": { + "additive_unit_conversion": false, + "bounds": [], + "data_path": "/Net/Groups/BGI/work_4/scratch/lalonso/GlobalForcingSet.zarr", + "depth_dimension": null, + "is_categorical": false, + "standard_name": null, + "sindbad_unit": null, + "source_product": "FLUXNET", + "source_to_sindbad_unit": 1.0, + "source_unit": null, + "source_variable": null, + "space_time_type": "spatiotemporal" + }, + "forcing_mask": { + "data_path": null, + "source_variable": null + }, + "variables": { + "f_ambient_CO2": { + "bounds": [200, 500], + "standard_name": "ambient_CO2", + "sindbad_unit": "ppm", + "source_unit": "ppm", + "source_variable": "atmCO2_SCRIPPS_global" + }, + "f_clay": { + "bounds": [0.0, 100.0], + "standard_name": "CLAY", + "sindbad_unit": "-", + "source_product": "soilgrids", + "source_to_sindbad_unit": 0.01, + "source_unit": "%", + "source_variable": "CLYPPT_soilgrid", + "space_time_type": "spatiovertical" + }, + "f_dist_intensity": { + "bounds": [-0.001, 1.1], + "is_categorical": true, + "standard_name": "isDisturbed flag for disturbance", + "sindbad_unit": null, + "source_product": "Besnard2018", + "source_unit": null, + "source_variable": "f_disturbance_zero2" + }, + "f_burnt_area": { + "bounds": [0.0, 1.0], + "standard_name": "fire fraction area per area pixel", + "sindbad_unit": null, + "source_unit": null, + "source_variable": "fire_frac_per_area" + }, + "f_tree_frac": { + "bounds": [0.0, 1.0], + "standard_name": "tree fraction", + "sindbad_unit": null, + "source_unit": null, + "source_variable": "tree_cover", + "space_time_type": "spatiovertical" + }, + "f_frac_vegetation": { + "bounds": [0.0, 1.0], + "standard_name": "tree fraction", + "sindbad_unit": null, + "source_unit": null, + "source_variable": "tree_cover", + "space_time_type": "spatiovertical" + }, + "f_pft": { + "bounds": [1, 17], + "standard_name": "pft index", + "sindbad_unit": null, + "source_unit": null, + "source_variable": "f_pft", + "space_time_type": "spatiovertical" + }, + "f_orgm": { + "bounds": [0.0, 100.0], + "standard_name": "CLAY", + "sindbad_unit": "-", + "source_product": "soilgrids", + "source_to_sindbad_unit": 0.0, + "source_unit": "%", + "source_variable": "OCSTHA_soilgrid", + "space_time_type": "spatiovertical" + }, + "f_PAR": { + "bounds": [0, 50], + "standard_name": "Photosynthetically active radiation (=Rg*.5)", + "sindbad_unit": "MJ m-2 d-1", + "source_to_sindbad_unit": 0.5, + "source_unit": "MJ m-2 d-1", + "source_variable": "ssrd" + }, + "f_rain": { + "bounds": [0, 300], + "standard_name": "Rain", + "sindbad_unit": "mm d-1", + "source_unit": "mm d-1", + "source_variable": "tp" + }, + "f_rg": { + "bounds": [0.0, 100.0], + "standard_name": "Global Radiation", + "sindbad_unit": "MJ m-2 d-1", + "source_unit": "MJ m-2 d-1", + "source_variable": "ssrd" + }, + "f_rg_pot": { + "bounds": [0.0, 100.0], + "standard_name": "Potential Global Radiation", + "sindbad_unit": "MJ m-2 d-1", + "source_unit": "MJ m-2 d-1", + "source_variable": "Rg_pot2" + }, + "f_rn": { + "bounds": [0.0, 100.0], + "standard_name": "Net radiation", + "sindbad_unit": "MJ m-2 d-1", + "source_unit": "MJ m-2 d-1", + "source_variable": "snr" + }, + "f_sand": { + "bounds": [0.0, 100.0], + "standard_name": "SAND", + "sindbad_unit": "-", + "source_product": "soilgrids", + "source_to_sindbad_unit": 0.01, + "source_unit": "%", + "source_variable": "SNDPPT_soilgrid", + "space_time_type": "spatiovertical" + }, + "f_silt": { + "bounds": [0.0, 100.0], + "standard_name": "CLAY", + "sindbad_unit": "-", + "source_product": "soilgrids", + "source_to_sindbad_unit": 0.01, + "source_unit": "%", + "source_variable": "SLTPPT_soilgrid", + "space_time_type": "spatiovertical" + }, + "f_airT": { + "additive_unit_conversion": true, + "bounds": [-80.0, 60.0], + "standard_name": "Tair", + "sindbad_unit": "°C", + "source_to_sindbad_unit": -273.15, + "source_unit": "K", + "source_variable": "t2m" + }, + "f_airT_day": { + "additive_unit_conversion": true, + "bounds": [-80.0, 60.0], + "standard_name": "TairDay", + "sindbad_unit": "°C", + "source_to_sindbad_unit": -273.15, + "source_unit": "K", + "source_variable": "t2m_daytime_mean" + }, + "f_VPD": { + "bounds": [0.01, 100.0], + "standard_name": "Vapor pressure deficit", + "sindbad_unit": "kPa", + "source_to_sindbad_unit": 0.1, + "source_unit": "hPa", + "source_variable": "vpd" + }, + "f_VPD_day": { + "bounds": [0.01, 100.0], + "standard_name": "Vapor pressure deficit Day", + "sindbad_unit": "kPa", + "source_to_sindbad_unit": 0.1, + "source_unit": "hPa", + "source_variable": "vpd_daytime_mean" + } + } + } \ No newline at end of file diff --git a/examples/exp_global_runs/settings_global_runs/forcing_f_backup.json b/examples/exp_global_runs/settings_global_runs/forcing_f_backup.json new file mode 100644 index 000000000..44a83397e --- /dev/null +++ b/examples/exp_global_runs/settings_global_runs/forcing_f_backup.json @@ -0,0 +1,184 @@ +{ + "data_dimension": { + "time": "Ti", + "permute": null, + "space": ["longitude", "latitude"] + }, + "default_forcing": { + "additive_unit_conversion": false, + "bounds": [], + "data_path": "/Net/Groups/BGI/work_4/scratch/lalonso/GlobalForcing.zarr", + "depth_dimension": null, + "is_categorical": false, + "standard_name": null, + "sindbad_unit": null, + "source_product": "FLUXNET", + "source_to_sindbad_unit": 1.0, + "source_unit": null, + "source_variable": null, + "space_time_type": "spatiotemporal" + }, + "forcing_mask": { + "data_path": null, + "source_variable": null + }, + "variables": { + "f_ambient_CO2": { + "bounds": [200, 500], + "standard_name": "ambient_CO2", + "sindbad_unit": "ppm", + "source_unit": "ppm", + "source_variable": "atmCO2_SCRIPPS_global" + }, + "f_clay": { + "bounds": [0.0, 100.0], + "standard_name": "CLAY", + "sindbad_unit": "-", + "source_product": "soilgrids", + "source_to_sindbad_unit": 0.01, + "source_unit": "%", + "source_variable": "CLYPPT_soilgrid", + "space_time_type": "spatiovertical" + }, + "f_dist_intensity": { + "bounds": [-0.001, 1.1], + "is_categorical": true, + "standard_name": "isDisturbed flag for disturbance", + "sindbad_unit": null, + "source_product": "Besnard2018", + "source_unit": null, + "source_variable": "f_disturbance_zero2" + }, + "f_burnt_area": { + "bounds": [0.0, 1.0], + "standard_name": "fire fraction area per area pixel", + "sindbad_unit": null, + "source_unit": null, + "source_variable": "fire_frac_per_area" + }, + "f_tree_frac": { + "bounds": [0.0, 1.0], + "standard_name": "tree fraction", + "sindbad_unit": null, + "source_unit": null, + "source_variable": "tree_cover", + "space_time_type": "spatiovertical" + }, + "f_frac_vegetation": { + "bounds": [0.0, 1.0], + "standard_name": "tree fraction", + "sindbad_unit": null, + "source_unit": null, + "source_variable": "tree_cover", + "space_time_type": "spatiovertical" + }, + "f_pft": { + "bounds": [1, 17], + "standard_name": "pft index", + "sindbad_unit": null, + "source_unit": null, + "source_variable": "f_pft", + "space_time_type": "spatiovertical" + }, + "f_orgm": { + "bounds": [0.0, 100.0], + "standard_name": "CLAY", + "sindbad_unit": "-", + "source_product": "soilgrids", + "source_to_sindbad_unit": 0.0, + "source_unit": "%", + "source_variable": "OCSTHA_soilgrid", + "space_time_type": "spatiovertical" + }, + "f_PAR": { + "bounds": [0, 50], + "standard_name": "Photosynthetically active radiation (=Rg*.5)", + "sindbad_unit": "MJ m-2 d-1", + "source_to_sindbad_unit": 0.5, + "source_unit": "MJ m-2 d-1", + "source_variable": "ssrd" + }, + "f_rain": { + "bounds": [0, 300], + "standard_name": "Rain", + "sindbad_unit": "mm d-1", + "source_unit": "mm d-1", + "source_variable": "tp" + }, + "f_rg": { + "bounds": [0.0, 100.0], + "standard_name": "Global Radiation", + "sindbad_unit": "MJ m-2 d-1", + "source_unit": "MJ m-2 d-1", + "source_variable": "ssrd" + }, + "f_rg_pot": { + "bounds": [0.0, 100.0], + "standard_name": "Potential Global Radiation", + "sindbad_unit": "MJ m-2 d-1", + "source_unit": "MJ m-2 d-1", + "source_variable": "Rg_pot2" + }, + "f_rn": { + "bounds": [0.0, 100.0], + "standard_name": "Net radiation", + "sindbad_unit": "MJ m-2 d-1", + "source_unit": "MJ m-2 d-1", + "source_variable": "snr" + }, + "f_sand": { + "bounds": [0.0, 100.0], + "standard_name": "SAND", + "sindbad_unit": "-", + "source_product": "soilgrids", + "source_to_sindbad_unit": 0.01, + "source_unit": "%", + "source_variable": "SNDPPT_soilgrid", + "space_time_type": "spatiovertical" + }, + "f_silt": { + "bounds": [0.0, 100.0], + "standard_name": "CLAY", + "sindbad_unit": "-", + "source_product": "soilgrids", + "source_to_sindbad_unit": 0.01, + "source_unit": "%", + "source_variable": "SLTPPT_soilgrid", + "space_time_type": "spatiovertical" + }, + "f_airT": { + "additive_unit_conversion": true, + "bounds": [-80.0, 60.0], + "standard_name": "Tair", + "sindbad_unit": "°C", + "source_to_sindbad_unit": -273.15, + "source_unit": "K", + "source_variable": "t2m" + }, + "f_airT_day": { + "additive_unit_conversion": true, + "bounds": [-80.0, 60.0], + "standard_name": "TairDay", + "sindbad_unit": "°C", + "source_to_sindbad_unit": -273.15, + "source_unit": "K", + "source_variable": "t2m_daytime_mean" + }, + "f_VPD": { + "bounds": [0.01, 100.0], + "standard_name": "Vapor pressure deficit", + "sindbad_unit": "kPa", + "source_to_sindbad_unit": 0.1, + "source_unit": "hPa", + "source_variable": "vpd" + }, + "f_VPD_day": { + "bounds": [0.01, 100.0], + "standard_name": "Vapor pressure deficit Day", + "sindbad_unit": "kPa", + "source_to_sindbad_unit": 0.1, + "source_unit": "hPa", + "source_variable": "vpd_daytime_mean" + } + } + } \ No newline at end of file diff --git a/examples/exp_global_runs/settings_global_runs/model_structure.json b/examples/exp_global_runs/settings_global_runs/model_structure.json new file mode 100644 index 000000000..01d647a62 --- /dev/null +++ b/examples/exp_global_runs/settings_global_runs/model_structure.json @@ -0,0 +1,258 @@ +{ + "default_model": { + "order": 0, + "use_in_spinup": true + }, + "models": { + "ambientCO2": { + "approach": "forcing" + }, + "autoRespiration": { + "approach": "Thornley2000A" + }, + "autoRespirationAirT": { + "approach": "Q10" + }, + "cAllocation": { + "approach": "GSI" + }, + "cAllocationLAI": { + "approach": "none" + }, + "cAllocationNutrients": { + "approach": "none" + }, + "cAllocationRadiation": { + "approach": "gpp" + }, + "cAllocationSoilT": { + "approach": "gpp" + }, + "cAllocationSoilW": { + "approach": "gpp" + }, + "cAllocationTreeFraction": { + "approach": "Friedlingstein1999" + }, + "capillaryFlow": { + "approach": "VanDijk2010" + }, + "cCycle": { + "approach": "GSI" + }, + "cCycleBase": { + "approach": "GSITOOPFT" + }, + "cCycleConsistency": { + "approach": "simple" + }, + "cVegetationDieOff": { + "approach": "forcing", + "use_in_spinup": false + }, + "cFireBurnedArea": { + "approach": "forcing", + "use_in_spinup": false + }, + "cFireCombustionCompleteness": { + "approach": "vanDerWerf2006", + "use_in_spinup": false + }, + "cFireMortality": { + "approach": "vanDerWerf2004", + "use_in_spinup": false + }, + "cCycleDisturbance": { + "approach": "WROASTEDTOO", + "use_in_spinup": false + }, + "cFlow": { + "approach": "GSI" + }, + "cFlowSoilProperties": { + "approach": "none" + }, + "cFlowVegProperties": { + "approach": "none" + }, + "cTau": { + "approach": "mult" + }, + "cTauLAI": { + "approach": "none" + }, + "cTauSoilProperties": { + "approach": "none" + }, + "cTauSoilT": { + "approach": "Q10" + }, + "cTauSoilW": { + "approach": "GSI" + }, + "cTauVegProperties": { + "approach": "none" + }, + "cBiomass": { + "approach": "treeGrass" + }, + "drainage": { + "approach": "dos" + }, + "evaporation": { + "approach": "fAPAR" + }, + "evapotranspiration": { + "approach": "sum" + }, + "vegFraction": { + "approach": "forcing" + }, + "fAPAR": { + "approach": "cVegLeafBareFrac" + }, + "getPools": { + "approach": "simple" + }, + "gpp": { + "approach": "coupled" + }, + "gppAirT": { + "approach": "CASA" + }, + "gppDemand": { + "approach": "mult" + }, + "gppDiffRadiation": { + "approach": "Wang2015" + }, + "gppDirRadiation": { + "approach": "none" + }, + "gppPotential": { + "approach": "Monteith" + }, + "gppSoilW": { + "approach": "Stocker2020" + }, + "gppVPD": { + "approach": "PRELES" + }, + "groundWRecharge": { + "approach": "dos" + }, + "groundWSoilWInteraction": { + "approach": "VanDijk2010" + }, + "LAI": { + "approach": "cVegLeaf" + }, + "percolation": { + "approach": "WBP" + }, + "PET": { + "approach": "Lu2005" + }, + "rainSnow": { + "approach": "Tair" + }, + "rootMaximumDepth": { + "approach": "fracSoilD" + }, + "rootWaterEfficiency": { + "approach": "expCvegRoot" + }, + "rootWaterUptake": { + "approach": "proportion" + }, + "runoff": { + "approach": "sum" + }, + "runoffBase": { + "approach": "Zhang2008" + }, + "runoffOverland": { + "approach": "Sat" + }, + "runoffSaturationExcess": { + "approach": "Bergstroem1992" + }, + "runoffSurface": { + "approach": "all" + }, + "snowFraction": { + "approach": "HTESSEL" + }, + "snowMelt": { + "approach": "TairRn" + }, + "soilProperties": { + "approach": "Saxton2006" + }, + "soilTexture": { + "approach": "forcing" + }, + "soilWBase": { + "approach": "uniform" + }, + "transpiration": { + "approach": "coupled" + }, + "transpirationSupply": { + "approach": "wAWC" + }, + "treeFraction": { + "approach": "forcing" + }, + "vegAvailableWater": { + "approach": "rootWaterEfficiency" + }, + "waterBalance": { + "approach": "simple", + "use_in_spinup": false + }, + "wCycle": { + "approach": "components" + }, + "wCycleBase": { + "approach": "simple" + }, + "WUE": { + "approach": "expVPDDayCo2" + } + }, + "pools": { + "carbon": { + "combine": "cEco", + "components": { + "cVeg": { + "Root": [1, 25.0], + "Wood": [1, 25.0], + "Leaf": [1, 25.0], + "Reserve": [1, 10.0] + }, + "cLit": { + "Fast": [1, 100.0], + "Slow": [1, 250.0] + }, + "cSoil": { + "Slow": [1, 500.0], + "Old": [1, 1000.0] + } + }, + "state_variables": {} + }, + "water": { + "combine": "TWS", + "components": { + "soilW": [[50.0, 200.0, 750.0, 1000.0], 100.0], + "groundW": [1, 1000.0], + "snowW": [1, 0.01], + "surfaceW": [1, 0.01] + }, + "state_variables": { + "Δ": 0.0 + } + } + } + } \ No newline at end of file From ca85f8999ee308ea0f358d50a46927be17545fa5 Mon Sep 17 00:00:00 2001 From: lazarusA Date: Mon, 21 Jul 2025 14:36:57 +0200 Subject: [PATCH 11/27] map params --- examples/exp_global_runs/Project.toml | 10 +++++ .../exp_global_runs/exp_global_parameters.jl | 38 +++++++++++++++++++ examples/exp_global_runs/exp_global_runs.jl | 13 ++++++- 3 files changed, 60 insertions(+), 1 deletion(-) create mode 100644 examples/exp_global_runs/Project.toml create mode 100644 examples/exp_global_runs/exp_global_parameters.jl diff --git a/examples/exp_global_runs/Project.toml b/examples/exp_global_runs/Project.toml new file mode 100644 index 000000000..34f6f8b80 --- /dev/null +++ b/examples/exp_global_runs/Project.toml @@ -0,0 +1,10 @@ +[deps] +Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" +Sindbad = "6686e6de-2010-4bf8-9178-b2bcc470766e" +SindbadData = "772b809f-5329-4bfd-a2e9-cc4714ce84b1" +SindbadML = "eb1d8004-2d42-4eef-b369-ed0645c101e8" +SindbadMetrics = "7a8077d0-ecf5-44b0-95d4-89eadbafc35c" +SindbadSetup = "2b7f2987-8a1e-48b4-8a02-cc9e7c4eeb1c" +SindbadTEM = "f6108451-10cb-42fa-b0c1-67671cf08f15" +SindbadUtils = "a2ad09c8-73f9-4387-ada0-5c5a8c58f8ab" +TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9" diff --git a/examples/exp_global_runs/exp_global_parameters.jl b/examples/exp_global_runs/exp_global_parameters.jl new file mode 100644 index 000000000..8247ca792 --- /dev/null +++ b/examples/exp_global_runs/exp_global_parameters.jl @@ -0,0 +1,38 @@ +ENV["JULIA_NUM_PRECOMPILE_TASKS"] = "1" +using SindbadTEM +using SindbadML +using SindbadData +using SindbadData.YAXArrays +using SindbadData.Zarr +using SindbadData.DimensionalData +using SindbadML.JLD2 +using TypedTables +using Flux + +# ? load trained neural network +# path_model = "/Net/Groups/BGI/work_4/scratch/lalonso/analysis/covs_all_fixed_rates/checkpoint/checkpoint_epoch_135.jld2" +path_model = joinpath(@__DIR__, "../analysis/modelFixK_ALL/checkpoint_epoch_500.jld2") + +trainedNN, lower_bound, upper_bound, ps_names, metadata_global = + loadTrainedNN(path_model) + +# ? load global PFT dataset +ds = open_dataset("/Net/Groups/BGI/work_5/scratch/lalonso/covariates_global_0d25.zarr") +cube_PFTs = readcubedata(ds["PFT_mask"]) +cube_KG = readcubedata(ds["KG_mask"]) + +ds_skip = open_dataset("/Net/Groups/BGI/work_5/scratch/lalonso/covariates_global_0d25.zarr"; + skip_keys=(:PFT_mask, :KG_mask),) + +ds_skip_cube = Cube(ds_skip) +c_keys = string.(ds_skip.cubes.keys) |> sort # ! this is a flaw, order should be consistent with training, hence it should be an output from there as well. + +ds_cubes_in = ds_skip_cube[Variables = At(c_keys)] +ds_cubes_in = readcubedata(ds_cubes_in) + +# ? compute and save new parameters +ps_path = "/Net/Groups/BGI/work_5/scratch/lalonso/parameters_ALL_0d25.zarr" + +out_params = mapParamsAll([cube_PFTs, cube_KG, ds_cubes_in], trainedNN, lower_bound, upper_bound, ps_names, ps_path; + metadata_global= metadata_global + ) \ No newline at end of file diff --git a/examples/exp_global_runs/exp_global_runs.jl b/examples/exp_global_runs/exp_global_runs.jl index 0c1faa85b..617d7225e 100644 --- a/examples/exp_global_runs/exp_global_runs.jl +++ b/examples/exp_global_runs/exp_global_runs.jl @@ -1 +1,12 @@ -2 + 2 \ No newline at end of file +ENV["JULIA_NUM_PRECOMPILE_TASKS"] = "1" +using SindbadTEM +using SindbadML +using SindbadData +using SindbadData.YAXArrays +using SindbadData.Zarr +using SindbadData.DimensionalData +using SindbadML.JLD2 +using TypedTables +using Flux + +## paths \ No newline at end of file From 3e2c37facb10a0825915f30a253aa55ecdca7d0b Mon Sep 17 00:00:00 2001 From: Lazaro Date: Mon, 21 Jul 2025 16:14:09 +0200 Subject: [PATCH 12/27] global parameters maps generation is working, tested only in Julia 1.11 --- examples/exp_global_runs/Project.toml | 2 ++ examples/exp_global_runs/exp_global_parameters.jl | 5 ++--- lib/SindbadData/Project.toml | 4 ++-- lib/SindbadML/Project.toml | 2 +- lib/SindbadML/src/SindbadML.jl | 2 +- 5 files changed, 8 insertions(+), 7 deletions(-) diff --git a/examples/exp_global_runs/Project.toml b/examples/exp_global_runs/Project.toml index 34f6f8b80..d585b020b 100644 --- a/examples/exp_global_runs/Project.toml +++ b/examples/exp_global_runs/Project.toml @@ -1,4 +1,5 @@ [deps] +DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" Sindbad = "6686e6de-2010-4bf8-9178-b2bcc470766e" SindbadData = "772b809f-5329-4bfd-a2e9-cc4714ce84b1" @@ -8,3 +9,4 @@ SindbadSetup = "2b7f2987-8a1e-48b4-8a02-cc9e7c4eeb1c" SindbadTEM = "f6108451-10cb-42fa-b0c1-67671cf08f15" SindbadUtils = "a2ad09c8-73f9-4387-ada0-5c5a8c58f8ab" TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9" +YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c" diff --git a/examples/exp_global_runs/exp_global_parameters.jl b/examples/exp_global_runs/exp_global_parameters.jl index 8247ca792..85d70b298 100644 --- a/examples/exp_global_runs/exp_global_parameters.jl +++ b/examples/exp_global_runs/exp_global_parameters.jl @@ -10,8 +10,7 @@ using TypedTables using Flux # ? load trained neural network -# path_model = "/Net/Groups/BGI/work_4/scratch/lalonso/analysis/covs_all_fixed_rates/checkpoint/checkpoint_epoch_135.jld2" -path_model = joinpath(@__DIR__, "../analysis/modelFixK_ALL/checkpoint_epoch_500.jld2") +path_model = joinpath("/Net/Groups/BGI/work_5/scratch/lalonso/checkpoint_epoch_208.jld2") trainedNN, lower_bound, upper_bound, ps_names, metadata_global = loadTrainedNN(path_model) @@ -31,7 +30,7 @@ ds_cubes_in = ds_skip_cube[Variables = At(c_keys)] ds_cubes_in = readcubedata(ds_cubes_in) # ? compute and save new parameters -ps_path = "/Net/Groups/BGI/work_5/scratch/lalonso/parameters_ALL_0d25.zarr" +ps_path = "/Net/Groups/BGI/work_5/scratch/lalonso/parameters_ALL_new_0d25.zarr" out_params = mapParamsAll([cube_PFTs, cube_KG, ds_cubes_in], trainedNN, lower_bound, upper_bound, ps_names, ps_path; metadata_global= metadata_global diff --git a/lib/SindbadData/Project.toml b/lib/SindbadData/Project.toml index 0478fdc35..2a5e8f7e9 100644 --- a/lib/SindbadData/Project.toml +++ b/lib/SindbadData/Project.toml @@ -25,6 +25,6 @@ DimensionalData = "0.29.12" FillArrays = "1.13.0" NCDatasets = "0.14.6" NetCDF = "0.12.1" -YAXArrayBase = "0.7.5" -YAXArrays = "0.6.0" +YAXArrayBase = "0.7" +YAXArrays = "0.6" Zarr = "0.9.4" diff --git a/lib/SindbadML/Project.toml b/lib/SindbadML/Project.toml index 5d2da08c0..d2448ae51 100644 --- a/lib/SindbadML/Project.toml +++ b/lib/SindbadML/Project.toml @@ -37,6 +37,6 @@ SindbadTEM = {path = "../../lib/SindbadTEM"} SindbadUtils = {path = "../../lib/SindbadUtils"} [compat] -DiskArrays = "0.4.11" +DiskArrays = "0.4" Enzyme = "0.13.38" PolyesterForwardDiff = "0.1.2" diff --git a/lib/SindbadML/src/SindbadML.jl b/lib/SindbadML/src/SindbadML.jl index ac6661b02..2cdd80168 100644 --- a/lib/SindbadML/src/SindbadML.jl +++ b/lib/SindbadML/src/SindbadML.jl @@ -86,5 +86,5 @@ module SindbadML include("siteLosses.jl") include("oneHots.jl") include("loadCovariates.jl") - + include("mapParams.jl") end From b9c60230922274bb916012e5be200eedd70e6081 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Tue, 22 Jul 2025 11:55:43 +0200 Subject: [PATCH 13/27] fire models by nuno --- examples/exp_global_runs/fire_models.jl | 86 +++++++++++++++++ src/Models/cFireBurnedArea/cFireBurnedArea.jl | 19 ++++ .../cFireBurnedArea_forcing.jl | 30 ++++++ .../cFireBurnedArea/cFireBurnedArea_none.jl | 12 +++ .../cFireCombustionCompleteness.jl | 20 ++++ .../cFireCombustionCompleteness_none.jl | 12 +++ ...reCombustionCompleteness_vanDerWerf2006.jl | 93 +++++++++++++++++++ src/Models/cFireMortality/cFireMortality.jl | 19 ++++ .../cFireMortality/cFireMortality_none.jl | 12 +++ .../cFireMortality_vanDerWerf2004.jl | 52 +++++++++++ 10 files changed, 355 insertions(+) create mode 100644 examples/exp_global_runs/fire_models.jl create mode 100644 src/Models/cFireBurnedArea/cFireBurnedArea.jl create mode 100644 src/Models/cFireBurnedArea/cFireBurnedArea_forcing.jl create mode 100644 src/Models/cFireBurnedArea/cFireBurnedArea_none.jl create mode 100644 src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness.jl create mode 100644 src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness_none.jl create mode 100644 src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness_vanDerWerf2006.jl create mode 100644 src/Models/cFireMortality/cFireMortality.jl create mode 100644 src/Models/cFireMortality/cFireMortality_none.jl create mode 100644 src/Models/cFireMortality/cFireMortality_vanDerWerf2004.jl diff --git a/examples/exp_global_runs/fire_models.jl b/examples/exp_global_runs/fire_models.jl new file mode 100644 index 000000000..413f4e625 --- /dev/null +++ b/examples/exp_global_runs/fire_models.jl @@ -0,0 +1,86 @@ +fire_models = (:constants, + :wCycleBase, + :rainSnow, + :rainIntensity, + :PET, + :ambientCO2, + :getPools, + :soilTexture, + :soilProperties, + :soilWBase, + :rootMaximumDepth, + :rootWaterEfficiency, + :PFT, + :plantForm, + :treeFraction, + :EVI, + :vegFraction, + :fAPAR, + :LAI, + :NDVI, + :NIRv, + :NDWI, + :snowFraction, + :sublimation, + :snowMelt, + :interception, + :runoffInfiltrationExcess, + :saturatedFraction, + :runoffSaturationExcess, + :runoffInterflow, + :runoffOverland, + :runoffSurface, + :runoffBase, + :percolation, + :evaporation, + :drainage, + :capillaryFlow, + :groundWRecharge, + :groundWSoilWInteraction, + :groundWSurfaceWInteraction, + :transpirationDemand, + :vegAvailableWater, + :transpirationSupply, + :gppPotential, + :gppDiffRadiation, + :gppDirRadiation, + :gppAirT, + :gppVPD, + :gppSoilW, + :gppDemand, + :WUE, + :gpp, + :transpiration, + :rootWaterUptake, + :cVegetationDieOff, + :cFireBurnedArea, + :cCycleBase, + :cCycleDisturbance, + :cTauSoilT, + :cTauSoilW, + :cTauLAI, + :cTauSoilProperties, + :cTauVegProperties, + :cTau, + :autoRespirationAirT, + :cAllocationLAI, + :cAllocationRadiation, + :cAllocationSoilW, + :cAllocationSoilT, + :cAllocationNutrients, + :cAllocation, + :cAllocationTreeFraction, + :autoRespiration, + :cFlowSoilProperties, + :cFlowVegProperties, + :cFlow, + :cCycleConsistency, + :cCycle, + :evapotranspiration, + :runoff, + :wCycle, + :waterBalance, + :cFireCombustionCompleteness, + :cFireMortality, + :cBiomass, + :deriveVariables) \ No newline at end of file diff --git a/src/Models/cFireBurnedArea/cFireBurnedArea.jl b/src/Models/cFireBurnedArea/cFireBurnedArea.jl new file mode 100644 index 000000000..b825c7c5f --- /dev/null +++ b/src/Models/cFireBurnedArea/cFireBurnedArea.jl @@ -0,0 +1,19 @@ +export cFireBurnedArea + +abstract type cFireBurnedArea <: LandEcosystem end + +include("cFireBurnedArea_none.jl") +include("cFireBurnedArea_forcing.jl") + +@doc """ +Accounts for carbon emissions due to fire + +# Approaches: +- none: no fire forcing, no emissions +- forcing: uses fire forcing data to calculate emissions + +*Created by* + - Nuno | nunocarvalhais + +""" +cFireBurnedArea \ No newline at end of file diff --git a/src/Models/cFireBurnedArea/cFireBurnedArea_forcing.jl b/src/Models/cFireBurnedArea/cFireBurnedArea_forcing.jl new file mode 100644 index 000000000..e85c167e4 --- /dev/null +++ b/src/Models/cFireBurnedArea/cFireBurnedArea_forcing.jl @@ -0,0 +1,30 @@ +export cFireBurnedArea_forcing + +struct cFireBurnedArea_forcing <: cFireBurnedArea end + +function define(params::cFireBurnedArea_forcing, forcing, land, helpers) + ## unpack land variables + @unpack_nt begin + f_burnt_area ⇐ forcing + end + c_fire_fba = f_burnt_area # in forcing called it f_burnt_area # f_burnt_area + # ## pack land variables + @pack_nt begin + c_fire_fba ⇒ land.diagnostics + end + return land +end + +function compute(params::cFireBurnedArea_forcing, forcing, land, helpers) + ## unpack land variables + @unpack_nt begin + f_burnt_area ⇐ forcing + end + # f_b = ... # some other pre process step, if needed + c_fire_fba = f_burnt_area # ? set to new value + + @pack_nt begin + c_fire_fba ⇒ land.diagnostics + end + return land +end \ No newline at end of file diff --git a/src/Models/cFireBurnedArea/cFireBurnedArea_none.jl b/src/Models/cFireBurnedArea/cFireBurnedArea_none.jl new file mode 100644 index 000000000..7a640f722 --- /dev/null +++ b/src/Models/cFireBurnedArea/cFireBurnedArea_none.jl @@ -0,0 +1,12 @@ +export cFireBurnedArea_none + +struct cFireBurnedArea_none <: cFireBurnedArea end + +function define(params::cFireBurnedArea_none, forcing, land, helpers) + return land +end + +@doc """ +TODO +""" +cFireBurnedArea_none diff --git a/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness.jl b/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness.jl new file mode 100644 index 000000000..83c817e12 --- /dev/null +++ b/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness.jl @@ -0,0 +1,20 @@ +export cFireCombustionCompleteness + +abstract type cFireCombustionCompleteness <: LandEcosystem end + +include("cFireCombustionCompleteness_vanDerWerf2006.jl") +include("cFireCombustionCompleteness_none.jl") + + +@doc """ +Accounts for carbon emissions due to fire, combustion completeness. + +# Approaches: +- none: no fire forcing, no emissions +- vanDerWerf2006: uses the van der Werf et al. (2006) method to calculate combustion completeness. + +*Created by* + - Nuno | nunocarvalhais + +""" +cFireCombustionCompleteness \ No newline at end of file diff --git a/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness_none.jl b/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness_none.jl new file mode 100644 index 000000000..8c4a79356 --- /dev/null +++ b/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness_none.jl @@ -0,0 +1,12 @@ +export cFireCombustionCompleteness_none + +struct cFireCombustionCompleteness_none <: cFireCombustionCompleteness end + +function define(params::cFireCombustionCompleteness_none, forcing, land, helpers) + return land +end + +@doc """ +TODO +""" +cFireCombustionCompleteness_none diff --git a/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness_vanDerWerf2006.jl b/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness_vanDerWerf2006.jl new file mode 100644 index 000000000..5d683b3a0 --- /dev/null +++ b/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness_vanDerWerf2006.jl @@ -0,0 +1,93 @@ +export cFireCombustionCompleteness_vanDerWerf2006 + +@with_kw struct cFireCombustionCompleteness_vanDerWerf2006{T1, T2, T3, T4, T5, T6, T7} <: cFireCombustionCompleteness + # fire combustion completeness parameters for stems, leaf; for the metabolic and structural parts of fine leaf litter (fcc_leaf_litter_m, fcc_leaf_litter_s)< for the coarse woody debris (literSlow) and for the organic soil layer (sol) + fcc_stem::T1 = [0.2f0, 0.3f0, 0.0f0, 1.0f0] # min, max, prev, current + fcc_leaf::T2 = [0.8f0, 1.0f0, 0.0f0, 1.0f0] + fcc_leaf_litter_m::T3 = [0.9f0, 1.0f0, 0.1f0, 0.9f0] + fcc_leaf_litter_s::T4 = [0.9f0, 1.0f0, 0.1f0, 0.9f0] + # fcc_sol::T5 = [0.9f0, 1.0f0, 0.1f0, 0.9f0] # we don't burnt organic soil + fcc_sol::T5 = [0.0f0, 0.0f0, 0.0f0, 1.0f0] + fcc_root::T6 = [0.0f0, 0.0f0, 0.0f0, 1.0f0] + fcc_cwd::T7 = [0.5f0, 0.6f0, 0.4f0, 0.6f0] +end + +function define(params::cFireCombustionCompleteness_vanDerWerf2006, forcing, land, helpers) + # @unpack_cFireCombustionCompleteness_vanDerWerf2006 params + ## instantiate variables + @unpack_nt begin + cEco ⇐ land.pools + zix ⇐ helpers.pools + end + + c_fire_ccMax = zero.(cEco) + c_fire_ccMin = zero.(cEco) + c_Fire_cci = zero.(cEco) + c_Fire_cc_fW = zero.(cEco) + # create CombustionCompletenessArray for each pool in zix + cc_lut = ( + :cVegRoot => :fcc_root, + :cVegWood => :fcc_stem, + :cVegReserve => :fcc_stem, + :cVegLeaf => :fcc_leaf, + :cLitFast => :fcc_leaf_litter_m, + :cLitSlow => :fcc_cwd, + :cSoilSlow => :fcc_sol, + ) + + for (k,v) in cc_lut + zix_keys = getproperty(zix, k) + imin, imax, _, _ = getproperty(params, v) # min, max, prev, current + # c_fire_ccMax[[zix_keys...]] .= imax + # c_fire_ccMin[[zix_keys...]] .= imin + for izix in zix_keys + @rep_elem imax ⇒ (c_fire_ccMax, izix, :cEco) + @rep_elem imin ⇒ (c_fire_ccMin, izix, :cEco) + end + end + + ## pack land variables + @pack_nt begin + (c_fire_ccMin, c_fire_ccMax, c_Fire_cci, c_Fire_cc_fW) ⇒ land.diagnostics + end + return land +end + +function compute(params::cFireCombustionCompleteness_vanDerWerf2006, forcing, land, helpers) + @unpack_cFireCombustionCompleteness_vanDerWerf2006 params + ## unpack land variables + @unpack_nt begin + (c_fire_ccMin, c_fire_ccMax, c_Fire_cci, c_Fire_cc_fW ) ⇐ land.diagnostics + gpp_f_soilW ⇐ land.diagnostics + zix ⇐ helpers.pools + soilW ⇐ land.pools + ∑w_sat ⇐ land.properties + (z_zero, o_one) ⇐ land.constants + end + + totalSoilW = maxZero(totalS(soilW)) + soilW_nor = minOne(totalSoilW / ∑w_sat) + + # for all soil pools c_Fire_cc_fW = soilW_nor + for zixSoil in (zix.cLit, zix.cSoil) + for izix in zixSoil + @rep_elem soilW_nor ⇒ (c_Fire_cc_fW, izix, :cEco) + end + end + # for all veg pools c_Fire_cc_fW = gpp_f_soilW + for zixVeg in zix.cVeg + @rep_elem gpp_f_soilW ⇒ (c_Fire_cc_fW, zixVeg, :cEco) + end + + # for all cEco pools + for zix_idx in zix.cEco + cci = (c_fire_ccMax[zix_idx] - c_fire_ccMin[zix_idx]) * (o_one - c_Fire_cc_fW[zix_idx]) + c_fire_ccMin[zix_idx] + @rep_elem cci ⇒ (c_Fire_cci, zix_idx, :cEco) + end + + # ## pack land variables + @pack_nt begin + (c_Fire_cci, c_Fire_cc_fW) ⇒ land.diagnostics + end + return land +end \ No newline at end of file diff --git a/src/Models/cFireMortality/cFireMortality.jl b/src/Models/cFireMortality/cFireMortality.jl new file mode 100644 index 000000000..081b429fb --- /dev/null +++ b/src/Models/cFireMortality/cFireMortality.jl @@ -0,0 +1,19 @@ +export cFireMortality + +abstract type cFireMortality <: LandEcosystem end + +include("cFireMortality_none.jl") +include("cFireMortality_vanDerWerf2004.jl") + + +@doc """ +Accounts for carbon emissions due to fire + +# Approaches: +- none +- vanDerWerf2004: uses the van der Werf et al. (2004) method to calculate fire mortality. + +*Created by* + - Nuno | nunocarvalhais +""" +cFireMortality \ No newline at end of file diff --git a/src/Models/cFireMortality/cFireMortality_none.jl b/src/Models/cFireMortality/cFireMortality_none.jl new file mode 100644 index 000000000..9432c9a64 --- /dev/null +++ b/src/Models/cFireMortality/cFireMortality_none.jl @@ -0,0 +1,12 @@ +export cFireMortality_none + +struct cFireMortality_none <: cFireMortality end + +function define(params::cFireMortality_none, forcing, land, helpers) + return land +end + +@doc """ +TODO +""" +cFireMortality_none diff --git a/src/Models/cFireMortality/cFireMortality_vanDerWerf2004.jl b/src/Models/cFireMortality/cFireMortality_vanDerWerf2004.jl new file mode 100644 index 000000000..1b214a04f --- /dev/null +++ b/src/Models/cFireMortality/cFireMortality_vanDerWerf2004.jl @@ -0,0 +1,52 @@ +export cFireMortality_vanDerWerf2004 + +@with_kw struct cFireMortality_vanDerWerf2004{T1, T2, T3, T4} <: cFireMortality + a::T1 = 0.01f0 + b::T2 = 0.59f0 + c::T3 = 0.6f0 + d::T4 = 0.25f0 +end + +function define(params::cFireMortality_vanDerWerf2004, forcing, land, helpers) + # @unpack_cFireMortality_vanDerWerf2006 params + ## instantiate variables + @unpack_nt begin + cEco ⇐ land.pools + end + c_Fire_k = one.(cEco) + ## pack land variables + @pack_nt begin + c_Fire_k ⇒ land.diagnostics + end + return land +end + +function compute(params::cFireMortality_vanDerWerf2004, forcing, land, helpers) + @unpack_cFireMortality_vanDerWerf2004 params + + @unpack_nt begin + c_Fire_k ⇐ land.diagnostics + frac_tree ⇐ land.states + zix ⇐ helpers.pools + (z_zero, o_one) ⇐ land.constants + end + # fire mortality according to Guido's paper + mortality = a + (b / (o_one + exp((c - frac_tree) * d))) + # wood mortality (wood / forest biomass lost) + for izix in zix.cVegWood + @rep_elem mortality ⇒ (c_Fire_k, izix, :cEco) + end + + # for the other vegetation pools the mortality scales with the frac_tree, we assume all the pools in grass have a mortality of 𝟙 + mortSplit = mortality * frac_tree + o_one * (o_one - frac_tree) + for c_izix in (zix.cVegRoot, zix.cVegLeaf, zix.cVegReserve) + for izix in c_izix + @rep_elem mortSplit ⇒ (c_Fire_k, izix, :cEco) + end + end + + @pack_nt begin + c_Fire_k ⇒ land.diagnostics + end + return land +end \ No newline at end of file From f899561af0c36c300a74a4848810d5480aee81d3 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Tue, 22 Jul 2025 15:41:31 +0200 Subject: [PATCH 14/27] follow the includeApproaches convention, but it would be better to simply use include --- examples/exp_global_runs/exp_global_runs.jl | 35 ++++++++++++---- src/Models/cFireBurnedArea/cFireBurnedArea.jl | 5 +-- .../cFireCombustionCompleteness.jl | 6 +-- src/Models/cFireMortality/cFireMortality.jl | 6 +-- src/Types/docStringForTypes.jl | 42 +++++++++++++++++++ 5 files changed, 75 insertions(+), 19 deletions(-) diff --git a/examples/exp_global_runs/exp_global_runs.jl b/examples/exp_global_runs/exp_global_runs.jl index 617d7225e..bdec7262e 100644 --- a/examples/exp_global_runs/exp_global_runs.jl +++ b/examples/exp_global_runs/exp_global_runs.jl @@ -1,12 +1,31 @@ -ENV["JULIA_NUM_PRECOMPILE_TASKS"] = "1" +# ENV["JULIA_NUM_PRECOMPILE_TASKS"] = "1" using SindbadTEM -using SindbadML using SindbadData -using SindbadData.YAXArrays -using SindbadData.Zarr -using SindbadData.DimensionalData +using SindbadSetup using SindbadML.JLD2 -using TypedTables -using Flux -## paths \ No newline at end of file +include("fire_models.jl"); + +experiment_output = "ALL_025" +remote_local = "/Net/Groups/BGI/tscratch/lalonso/" +remote_root = joinpath(remote_local, experiment_output) +mkpath(remote_root) + +domain = "Global"; +replace_info_spatial = Dict( + "experiment.model_output.path" => remote_root, + "experiment.basics.domain" => domain * "_spatial", + "experiment.flags.spinup_TEM" => true, + "experiment.flags.run_optimization" => false, + "info.settings.experiment.flags.calc_cost" => false, + "model_structure.sindbad_models" => fire_models + ); + +experiment_json = "../exp_global_runs/settings_global_runs/experiment.json"; + +info = getExperimentInfo(experiment_json; replace_info=replace_info_spatial); + +path_model = joinpath("/Net/Groups/BGI/work_5/scratch/lalonso/checkpoint_epoch_208.jld2") +model_props = JLD2.load(path_model) +tbl_params = model_props["parameter_table"] + diff --git a/src/Models/cFireBurnedArea/cFireBurnedArea.jl b/src/Models/cFireBurnedArea/cFireBurnedArea.jl index b825c7c5f..90a465f9a 100644 --- a/src/Models/cFireBurnedArea/cFireBurnedArea.jl +++ b/src/Models/cFireBurnedArea/cFireBurnedArea.jl @@ -1,9 +1,8 @@ export cFireBurnedArea abstract type cFireBurnedArea <: LandEcosystem end - -include("cFireBurnedArea_none.jl") -include("cFireBurnedArea_forcing.jl") +purpose(::Type{cFireBurnedArea}) = "Disturbance of the carbon cycle pools due to fire." +includeApproaches(cFireBurnedArea, @__DIR__) @doc """ Accounts for carbon emissions due to fire diff --git a/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness.jl b/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness.jl index 83c817e12..558484cc8 100644 --- a/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness.jl +++ b/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness.jl @@ -1,10 +1,8 @@ export cFireCombustionCompleteness abstract type cFireCombustionCompleteness <: LandEcosystem end - -include("cFireCombustionCompleteness_vanDerWerf2006.jl") -include("cFireCombustionCompleteness_none.jl") - +purpose(::Type{cFireCombustionCompleteness}) = "Disturbance of the carbon cycle pools due to fire, combustion completeness." +includeApproaches(cFireCombustionCompleteness, @__DIR__) @doc """ Accounts for carbon emissions due to fire, combustion completeness. diff --git a/src/Models/cFireMortality/cFireMortality.jl b/src/Models/cFireMortality/cFireMortality.jl index 081b429fb..73245451f 100644 --- a/src/Models/cFireMortality/cFireMortality.jl +++ b/src/Models/cFireMortality/cFireMortality.jl @@ -1,10 +1,8 @@ export cFireMortality abstract type cFireMortality <: LandEcosystem end - -include("cFireMortality_none.jl") -include("cFireMortality_vanDerWerf2004.jl") - +purpose(::Type{cFireMortality}) = "Disturbance of the carbon cycle pools due to fire." +includeApproaches(cFireMortality, @__DIR__) @doc """ Accounts for carbon emissions due to fire diff --git a/src/Types/docStringForTypes.jl b/src/Types/docStringForTypes.jl index c09fed978..9670a57d7 100644 --- a/src/Types/docStringForTypes.jl +++ b/src/Types/docStringForTypes.jl @@ -3831,3 +3831,45 @@ Use Optimisers.jl Descent optimizer for training ML models in SINDBAD """ Sindbad.Types.OptimisersDescent +@doc """ + +# cFireBurnedArea_forcing + +nothing + +## Type Hierarchy + +```cFireBurnedArea_forcing <: cFireBurnedArea <: LandEcosystem <: ModelTypes <: SindbadTypes <: Any``` + + +""" +Sindbad.Models.cFireBurnedArea_forcing + +@doc """ + +# cFireCombustionCompleteness_vanDerWerf2006 + +nothing + +## Type Hierarchy + +```cFireCombustionCompleteness_vanDerWerf2006 <: cFireCombustionCompleteness <: LandEcosystem <: ModelTypes <: SindbadTypes <: Any``` + + +""" +Sindbad.Models.cFireCombustionCompleteness_vanDerWerf2006 + +@doc """ + +# cFireMortality_vanDerWerf2004 + +nothing + +## Type Hierarchy + +```cFireMortality_vanDerWerf2004 <: cFireMortality <: LandEcosystem <: ModelTypes <: SindbadTypes <: Any``` + + +""" +Sindbad.Models.cFireMortality_vanDerWerf2004 + From 4c47cc26a40ab08a33b56caacbe77c890c3923ad Mon Sep 17 00:00:00 2001 From: Lazaro Date: Tue, 22 Jul 2025 16:24:19 +0200 Subject: [PATCH 15/27] use plantForm, also a bug, experiment should run also without optimization file, which is not the case --- .../settings_global_runs/experiment.json | 3 +- .../settings_global_runs/model_structure.json | 6 +- .../settings_global_runs/optimization.json | 258 ++++++++++++++++++ 3 files changed, 263 insertions(+), 4 deletions(-) create mode 100644 examples/exp_global_runs/settings_global_runs/optimization.json diff --git a/examples/exp_global_runs/settings_global_runs/experiment.json b/examples/exp_global_runs/settings_global_runs/experiment.json index 07e9880b4..971bfc0db 100644 --- a/examples/exp_global_runs/settings_global_runs/experiment.json +++ b/examples/exp_global_runs/settings_global_runs/experiment.json @@ -2,7 +2,8 @@ "basics": { "config_files": { "forcing": "forcing.json", - "model_structure": "model_structure.json" + "model_structure": "model_structure.json", + "optimization": "optimization.json" }, "domain": "Global", "name": "Fix_ALL_ALL", diff --git a/examples/exp_global_runs/settings_global_runs/model_structure.json b/examples/exp_global_runs/settings_global_runs/model_structure.json index 01d647a62..b903b17f4 100644 --- a/examples/exp_global_runs/settings_global_runs/model_structure.json +++ b/examples/exp_global_runs/settings_global_runs/model_structure.json @@ -1,6 +1,6 @@ { "default_model": { - "order": 0, + "implicit_t_repeat": 1, "use_in_spinup": true }, "models": { @@ -41,7 +41,7 @@ "approach": "GSI" }, "cCycleBase": { - "approach": "GSITOOPFT" + "approach": "GSI_PlantForm" }, "cCycleConsistency": { "approach": "simple" @@ -63,7 +63,7 @@ "use_in_spinup": false }, "cCycleDisturbance": { - "approach": "WROASTEDTOO", + "approach": "WROASTED", "use_in_spinup": false }, "cFlow": { diff --git a/examples/exp_global_runs/settings_global_runs/optimization.json b/examples/exp_global_runs/settings_global_runs/optimization.json new file mode 100644 index 000000000..36271db8e --- /dev/null +++ b/examples/exp_global_runs/settings_global_runs/optimization.json @@ -0,0 +1,258 @@ +{ + "algorithm": "CMAEvolutionStrategy_CMAES", + "model_parameter_default": { + "distribution": ["normal", [0.0, 1.0]], + "is_ml": true + }, + "model_parameters_to_optimize": { + "autoRespiration,RMN": null, + "autoRespiration,YG": null, + "autoRespirationAirT,Q10": null, + "cCycleBase,c_remain": null, + "cCycleBase,ηH": null, + "cFlow,f_τ": null, + "cFlow,k_shedding": null, + "cFlow,slope_leaf_root_to_reserve": null, + "cFlow,slope_reserve_to_leaf_root": null, + "cTauSoilT,Q10": null, + "cTauSoilW,opt_soilW": null, + "drainage,dos_exp": null, + "evaporation,k_evaporation": null, + "evaporation,α": null, + "fAPAR,k_extinction": null, + "gppAirT,opt_airT": null, + "gppAirT,opt_airT_A": null, + "gppAirT,opt_airT_B": null, + "gppDiffRadiation,μ": null, + "gppPotential,εmax": null, + "gppSoilW,q": null, + "gppSoilW,θstar": null, + "gppVPD,sat_ambient_CO2": null, + "gppVPD,c_κ": null, + "gppVPD,κ": null, + "groundWRecharge,dos_exp": null, + "rootMaximumDepth,constant_frac_max_root_depth": null, + "rootWaterEfficiency,k_efficiency_cVegRoot": null, + "runoffSaturationExcess,β": null, + "snowMelt,melt_Rn": null, + "snowMelt,melt_T": null, + "transpirationSupply,k_transpiration": null, + "WUE,WUE_one_hpa": null, + "WUE,κ": null + }, + "multi_constraint_method": "cost_sum", + "multi_objective_algorithm": false, + "observational_constraints": [ + "gpp", + "nee", + "reco", + "transpiration", + "evapotranspiration", + "agb", + "ndvi" + ], + "observations": { + "default_cost": { + "aggr_func": "nanmean", + "aggr_obs": false, + "aggr_order": "time_space", + "cost_metric": "NNSE_inv", + "cost_weight": 1.0, + "min_data_points": 1, + "spatial_data_aggr": "concat_data", + "spatial_cost_aggr": "spatially_variable", + "spatial_weight": false, + "temporal_data_aggr": "day" + }, + "default_observation": { + "additive_unit_conversion": false, + "bounds": null, + "data_path": "../data/FLUXNET_v2023_12_1D.zarr", + "is_categorical": false, + "standard_name": null, + "sindbad_unit": null, + "source_product": "", + "source_to_sindbad_unit": 1.0, + "source_unit": null, + "source_variable": null, + "space_time_type": "spatiotemporal" + }, + "use_quality_flag": true, + "use_spatial_weight": false, + "use_uncertainty": false, + "variables": { + "agb": { + "cost_options": { + "cost_metric": "NMAE1R" + }, + "data": { + "bounds": [0.0, 100000.0], + "standard_name": "Above-ground biomass", + "sindbad_unit": "gC m-2", + "source_unit": "gC m-2", + "source_variable": "agb_merged_PFT" + }, + "model_full_var": "states.aboveground_biomass" + }, + "evapotranspiration": { + "cost_options": { + "cost_metric": "NNSE_inv" + }, + "data": { + "bounds": [0.0, 100.0], + "is_categorical": false, + "standard_name": "Evapotranspiration", + "sindbad_unit": "mm day-1", + "source_to_sindbad_unit": 0.4081632653, + "source_unit": "MJ m-2 d-1", + "source_variable": "LE" + }, + "model_full_var": "fluxes.evapotranspiration", + "qflag": { + "bounds": [0.85, 1.0], + "data_path": null, + "source_variable": "LE_QC_merged" + }, + "unc": { + "bounds": [0.0, 100.0], + "data_path": null, + "source_to_sindbad_unit": 0.4081632653, + "source_variable": "LE_RANDUNC" + } + }, + "gpp": { + "cost_options": { + "cost_metric": "NNSE_inv" + }, + "data": { + "bounds": [0.0, 100.0], + "standard_name": "Gross Primary Productivity", + "sindbad_unit": "gC m-2 day-1", + "source_unit": "gC m-2 day-1", + "source_variable": "GPP_NT" + }, + "model_full_var": "fluxes.gpp", + "qflag": { + "bounds": [0.85, 1.0], + "data_path": null, + "source_variable": "GPP_QC_NT_merged" + }, + "unc": { + "bounds": [0.0, 400.0], + "data_path": null, + "source_variable": "NEE_RANDUNC" + } + }, + "ndvi": { + "cost_options": { + "cost_metric": "NNSE_inv", + "temporal_data_aggr": "day_anomaly", + "aggr_obs": true + }, + "data": { + "bounds": [0.0, 1.0], + "standard_name": "NDVI_MCD43A", + "sindbad_unit": "adimensional", + "source_unit": "adimensional", + "source_variable": "NDVI_MCD43A" + }, + "model_full_var": "states.fAPAR_bare", + "qflag": { + "bounds": [0.85, 1.0], + "data_path": null, + "source_variable": "NDVI_QC_MCD43A" + } + }, + "nee": { + "cost_options": { + "cost_metric": "NNSE_inv" + }, + "data": { + "bounds": [-100.0, 100.0], + "standard_name": "Net Ecosystem Exchange", + "sindbad_unit": "gC m-2 day-1", + "source_unit": "gC m-2 day-1", + "source_variable": "NEE" + }, + "model_full_var": "fluxes.nee", + "qflag": { + "bounds": [0.85, 1.0], + "data_path": null, + "source_variable": "NEE_QC_merged" + }, + "unc": { + "bounds": [0.0, 400.0], + "data_path": null, + "source_variable": "NEE_RANDUNC" + } + }, + "nirv": { + "cost_options": { + "cost_metric": "n_pcor_inv", + "temporal_data_aggr": "day_anomaly", + "aggr_obs": true + }, + "data": { + "bounds": [0.0, 1.0], + "standard_name": "Near Infrared Reflectance", + "sindbad_unit": "adimensional", + "source_unit": "-", + "source_variable": "NIRv_MCD43A" + }, + "model_full_var": "fluxes.gpp", + "qflag": { + "bounds": [0.85, 1.0], + "data_path": null, + "source_variable": "NIRv_QC_MCD43A" + } + }, + "reco": { + "cost_options": { + "cost_metric": "NNSE_inv" + }, + "data": { + "bounds": [0.0, 100.0], + "standard_name": "Ecosystem Respiration", + "sindbad_unit": "gC m-2 day-1", + "source_unit": "gC m-2 day-1", + "source_variable": "RECO_NT" + }, + "model_full_var": "fluxes.eco_respiration", + "qflag": { + "bounds": [0.85, 1.0], + "data_path": null, + "source_variable": "RECO_QC_NT_merged" + }, + "unc": { + "bounds": [0.0, 400.0], + "data_path": null, + "source_variable": "NEE_RANDUNC" + } + }, + "transpiration": { + "cost_options": { + "cost_metric": "NNSE_inv" + }, + "data": { + "bounds": [0.0, 100.0], + "standard_name": "transpiration", + "sindbad_unit": "mm day-1", + "source_unit": "mm day-1", + "source_variable": "T_NT_TEA" + }, + "model_full_var": "fluxes.transpiration", + "qflag": { + "bounds": [0.85, 1.0], + "data_path": null, + "source_variable": "LE_QC_merged" + }, + "unc": { + "bounds": [0.0, 100.0], + "data_path": null, + "source_to_sindbad_unit": 0.4081632653, + "source_variable": "LE_RANDUNC" + } + } + } + } +} \ No newline at end of file From 64d5144464abd1e96d41ab7b81a458bfc42e9619 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Thu, 24 Jul 2025 11:42:17 +0200 Subject: [PATCH 16/27] fixes getForcing dispatch for YAXArray incubes --- examples/exp_global_runs/Project.toml | 1 + examples/exp_global_runs/exp_global_runs.jl | 4 +++ lib/SindbadData/src/SindbadData.jl | 1 + lib/SindbadData/src/getForcing.jl | 31 +++++++++++++++++---- lib/SindbadData/src/utilsData.jl | 2 +- 5 files changed, 33 insertions(+), 6 deletions(-) diff --git a/examples/exp_global_runs/Project.toml b/examples/exp_global_runs/Project.toml index d585b020b..59ae968f0 100644 --- a/examples/exp_global_runs/Project.toml +++ b/examples/exp_global_runs/Project.toml @@ -10,3 +10,4 @@ SindbadTEM = "f6108451-10cb-42fa-b0c1-67671cf08f15" SindbadUtils = "a2ad09c8-73f9-4387-ada0-5c5a8c58f8ab" TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9" YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c" +Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99" diff --git a/examples/exp_global_runs/exp_global_runs.jl b/examples/exp_global_runs/exp_global_runs.jl index bdec7262e..fa3d78e41 100644 --- a/examples/exp_global_runs/exp_global_runs.jl +++ b/examples/exp_global_runs/exp_global_runs.jl @@ -29,3 +29,7 @@ path_model = joinpath("/Net/Groups/BGI/work_5/scratch/lalonso/checkpoint_epoch_2 model_props = JLD2.load(path_model) tbl_params = model_props["parameter_table"] +forcing = getForcing(info); +ds = open_dataset("/Net/Groups/BGI/work_4/scratch/lalonso/GlobalForcingSet.zarr"); + +run_helpers = prepTEM(forcing, info); \ No newline at end of file diff --git a/lib/SindbadData/src/SindbadData.jl b/lib/SindbadData/src/SindbadData.jl index faac7cc24..d745df9dc 100644 --- a/lib/SindbadData/src/SindbadData.jl +++ b/lib/SindbadData/src/SindbadData.jl @@ -43,6 +43,7 @@ module SindbadData using AxisKeys: KeyedArray, AxisKeys using FillArrays using DimensionalData + using DimensionalData: DimensionalData as DD using NCDatasets @reexport using NetCDF @reexport using YAXArrays diff --git a/lib/SindbadData/src/getForcing.jl b/lib/SindbadData/src/getForcing.jl index 0e695f3c4..6b1fed3b9 100644 --- a/lib/SindbadData/src/getForcing.jl +++ b/lib/SindbadData/src/getForcing.jl @@ -97,11 +97,7 @@ function createForcingNamedTuple(incubes, f_sizes, f_dimensions, info) typed_cubes = getInputArrayOfType(incubes, input_array_type) data_ts_type=[] for incube in typed_cubes - if in(:time, AxisKeys.dimnames(incube)) - push!(data_ts_type, ForcingWithTime()) - else - push!(data_ts_type, ForcingWithoutTime()) - end + push!(data_ts_type, collectTypesTiDim(incube)) end data_ts_type = [_dt for _dt in data_ts_type] f_types = Tuple(Tuple.(Pair.(forcing_vars, data_ts_type))) @@ -115,6 +111,30 @@ function createForcingNamedTuple(incubes, f_sizes, f_dimensions, info) return forcing end +# Fixes bug. Collecting forcing variables `WithTime` was not working for YAXArrays +# These two functions now dispatch on KeyedArray and YAXArray +""" + collectTypesTiDim(incube::AxisKeys.KeyedArray) +""" +function collectTypesTiDim(incube::AxisKeys.KeyedArray) + if in(:time, AxisKeys.dimnames(incube)) + return ForcingWithTime() + else + return ForcingWithoutTime() + end +end + +""" + collectTypesTiDim(incube::YAXArrays.YAXArray) +""" +function collectTypesTiDim(incube::YAXArrays.YAXArray) + if hasdim(incube, Ti) || hasdim(incube, DD.Dim{:Time}) || hasdim(incube, DD.Dim{:time}) + return ForcingWithTime() + else + return ForcingWithoutTime() + end +end + """ getForcing(info::NamedTuple) @@ -160,6 +180,7 @@ function getForcing(info::NamedTuple) default_info = info.experiment.data_settings.forcing.default_forcing forcing_vars = keys(forcing_data_settings.variables) + @show forcing_vars tar_dims = getTargetDimensionOrder(info) showInfo(getForcing, @__FILE__, @__LINE__, "getting forcing variables. Units given in forcing settings are not strictly enforced but shown for reference. Bounds are applied after unit conversion...", n_m=1) vinfo = nothing diff --git a/lib/SindbadData/src/utilsData.jl b/lib/SindbadData/src/utilsData.jl index 5ea778a8c..bab729110 100644 --- a/lib/SindbadData/src/utilsData.jl +++ b/lib/SindbadData/src/utilsData.jl @@ -430,7 +430,7 @@ function subsetAndProcessYax(yax, forcing_mask, tar_dims, _data_info, info, ::Va end #todo mean of the data instead of zero or nan - if info.settings.experiment.model_output.output_array_type == "array" + if info.experiment.exe_rules.land_output_type == "array" vfill = 0.0 if fill_nan vfill = NaN From e55c916b175b6c64fd374f51f279bc82482a735a Mon Sep 17 00:00:00 2001 From: Lazaro Date: Thu, 24 Jul 2025 14:09:51 +0200 Subject: [PATCH 17/27] :'( more updates to internals, now global forcings with lazy loading works --- examples/exp_global_runs/exp_global_runs.jl | 3 +-- lib/SindbadData/src/getForcing.jl | 27 +++++++++++++++++---- lib/SindbadTEM/src/prepTEMOut.jl | 2 +- 3 files changed, 24 insertions(+), 8 deletions(-) diff --git a/examples/exp_global_runs/exp_global_runs.jl b/examples/exp_global_runs/exp_global_runs.jl index fa3d78e41..0150f99cd 100644 --- a/examples/exp_global_runs/exp_global_runs.jl +++ b/examples/exp_global_runs/exp_global_runs.jl @@ -30,6 +30,5 @@ model_props = JLD2.load(path_model) tbl_params = model_props["parameter_table"] forcing = getForcing(info); -ds = open_dataset("/Net/Groups/BGI/work_4/scratch/lalonso/GlobalForcingSet.zarr"); - +# ds = open_dataset("/Net/Groups/BGI/work_4/scratch/lalonso/GlobalForcingSet.zarr"); run_helpers = prepTEM(forcing, info); \ No newline at end of file diff --git a/lib/SindbadData/src/getForcing.jl b/lib/SindbadData/src/getForcing.jl index 6b1fed3b9..0cf4be475 100644 --- a/lib/SindbadData/src/getForcing.jl +++ b/lib/SindbadData/src/getForcing.jl @@ -21,11 +21,7 @@ function collectForcingSizes(info, in_yax) dnames = Symbol[] dsizes = [] push!(dnames, time_dim_name) - if time_dim_name in in_yax - push!(dsizes, length(getproperty(in_yax, time_dim_name))) - else - push!(dsizes, length(DimensionalData.lookup(in_yax, time_dim_name))) - end + collectTimeSizes!(in_yax, dsizes, time_dim_name) for space ∈ info.experiment.data_settings.forcing.data_dimension.space push!(dnames, Symbol(space)) push!(dsizes, length(getproperty(in_yax, Symbol(space)))) @@ -34,6 +30,27 @@ function collectForcingSizes(info, in_yax) return f_sizes end + +function collectTimeSizes!(incube::YAXArrays.YAXArray, dsizes, time_dim_name) + new_size = length(DimensionalData.lookup(incube, time_dim_name)) + push!(dsizes, new_size) + return nothing +end + +function collectTimeSizes!(incube::AxisKeys.KeyedArray, dsizes, time_dim_name) + new_size = length(getproperty(incube, time_dim_name)) + push!(dsizes, new_size) + return nothing +end + +""" + collectTimeSizes!(incube::YAXArrays.YAXArray, dsizes, time_dim_name) + collectTimeSizes!(incube::AxisKeys.KeyedArray, dsizes, time_dim_name) + +Collects the size of the time dimension from the input cube and appends it to `dsizes`. +""" +function collectTimeSizes! end + """ collectForcingHelpers(info, f_sizes, f_dimensions) diff --git a/lib/SindbadTEM/src/prepTEMOut.jl b/lib/SindbadTEM/src/prepTEMOut.jl index 679e63d12..34efe3d13 100644 --- a/lib/SindbadTEM/src/prepTEMOut.jl +++ b/lib/SindbadTEM/src/prepTEMOut.jl @@ -160,7 +160,7 @@ function getOutDims(info, forcing_helpers, ::OutputYAXArray) vname = string(last(vname_full)) _properties = collectMetadata(info, vname_full) vdims = var_dims[v_index] - outformat = info.settings.experiment.model_output.format + outformat = info.output.format backend = outformat == "nc" ? :netcdf : :zarr out_dim = YAXArrays.OutDims(vdims...; properties = _properties, From 1cabf9f13e43785e7848afa265e0bfad68b99c64 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Thu, 24 Jul 2025 16:53:09 +0200 Subject: [PATCH 18/27] missing exports, internals and rm application of getSpinupTemLite, it should come last --- examples/exp_global_runs/exp_global_runs.jl | 22 ++++++++++++++++++++- lib/SindbadData/src/utilsData.jl | 3 +++ lib/SindbadTEM/src/SindbadTEM.jl | 5 ++++- lib/SindbadTEM/src/prepTEM.jl | 4 +++- lib/SindbadTEM/src/runTEMCube.jl | 10 +++++----- 5 files changed, 36 insertions(+), 8 deletions(-) diff --git a/examples/exp_global_runs/exp_global_runs.jl b/examples/exp_global_runs/exp_global_runs.jl index 0150f99cd..dc38241ab 100644 --- a/examples/exp_global_runs/exp_global_runs.jl +++ b/examples/exp_global_runs/exp_global_runs.jl @@ -1,4 +1,6 @@ # ENV["JULIA_NUM_PRECOMPILE_TASKS"] = "1" +using Sindbad +using SindbadUtils using SindbadTEM using SindbadData using SindbadSetup @@ -31,4 +33,22 @@ tbl_params = model_props["parameter_table"] forcing = getForcing(info); # ds = open_dataset("/Net/Groups/BGI/work_4/scratch/lalonso/GlobalForcingSet.zarr"); -run_helpers = prepTEM(forcing, info); \ No newline at end of file +run_helpers = prepTEM(forcing, info); + +tuple_models = getTupleFromLongTuple(info.models.forward); + +# newly generated parameters +ps_path = "/Net/Groups/BGI/work_5/scratch/lalonso/parameters_ALL_new_0d25.zarr" + +in_cube_params = Cube(ps_path) +in_cube_ps = permutedims(in_cube_params, (2,3,1)) +in_cube_ps = readcubedata(in_cube_ps) + +yax_max_cache = info.experiment.exe_rules.yax_max_cache + +outcubes = SindbadTEM.runTEMYaxParameters( + tuple_models, + forcing, + in_cube_ps, + tbl_params, + info); \ No newline at end of file diff --git a/lib/SindbadData/src/utilsData.jl b/lib/SindbadData/src/utilsData.jl index bab729110..636d0005e 100644 --- a/lib/SindbadData/src/utilsData.jl +++ b/lib/SindbadData/src/utilsData.jl @@ -1,5 +1,8 @@ export getNumberOfTimeSteps export mapCleanData +export cleanData +export applyUnitConversion +export applyQCBound export subsetAndProcessYax export yaxCubeToKeyedArray export toDimStackArray diff --git a/lib/SindbadTEM/src/SindbadTEM.jl b/lib/SindbadTEM/src/SindbadTEM.jl index 101eed0ae..e3332a329 100644 --- a/lib/SindbadTEM/src/SindbadTEM.jl +++ b/lib/SindbadTEM/src/SindbadTEM.jl @@ -58,9 +58,12 @@ module SindbadTEM using Sindbad using SindbadUtils using SindbadSetup - using SindbadData: YAXArrays using ThreadPools + using SindbadData.YAXArrays + using SindbadData.Zarr + using SindbadData: AllNaN + using SindbadData: cleanData include("utilsTEM.jl") include("deriveSpinupForcing.jl") diff --git a/lib/SindbadTEM/src/prepTEM.jl b/lib/SindbadTEM/src/prepTEM.jl index 184e6bb33..44b54a2bd 100644 --- a/lib/SindbadTEM/src/prepTEM.jl +++ b/lib/SindbadTEM/src/prepTEM.jl @@ -190,7 +190,9 @@ function getRunTEMInfo(info, forcing) upd_tem_helpers = setTupleField(upd_tem_helpers, (:vals, vals)) upd_tem_helpers = setTupleField(upd_tem_helpers, (:model_helpers, model_helpers)) upd_tem_helpers = setTupleField(upd_tem_helpers, (:run, tem_helpers.run)) - upd_tem_helpers = setTupleField(upd_tem_helpers, (:spinup_sequence, getSpinupTemLite(info.spinup.sequence))) + # upd_tem_helpers = setTupleField(upd_tem_helpers, (:spinup_sequence, getSpinupTemLite(info.spinup.sequence))) + # don't do `getSpinupTemLite` here, but rather later on more inner functions! + upd_tem_helpers = setTupleField(upd_tem_helpers, (:spinup_sequence, info.spinup.sequence)) return upd_tem_helpers end diff --git a/lib/SindbadTEM/src/runTEMCube.jl b/lib/SindbadTEM/src/runTEMCube.jl index db72fd4e5..73b141581 100644 --- a/lib/SindbadTEM/src/runTEMCube.jl +++ b/lib/SindbadTEM/src/runTEMCube.jl @@ -135,10 +135,10 @@ end """ - run_psTEMYax(selected_models::Tuple, forcing::NamedTuple, in_cube_params, tbl_params, info::NamedTuple) + runTEMYaxParameters(selected_models::Tuple, forcing::NamedTuple, in_cube_params, tbl_params, info::NamedTuple) """ -function run_psTEMYax(selected_models::Tuple, forcing::NamedTuple, in_cube_params, tbl_params, info::NamedTuple) +function runTEMYaxParameters(selected_models::Tuple, forcing::NamedTuple, in_cube_params, tbl_params, info::NamedTuple) # forcing/input information in_cubes_forcing = forcing.data @@ -151,9 +151,9 @@ function run_psTEMYax(selected_models::Tuple, forcing::NamedTuple, in_cube_param loc_land = deepcopy(run_helpers.loc_land) _data_fill = 0.0f0 - _forcing_default_info = info.settings.forcing.default_forcing + _forcing_default_info = info.experiment.data_settings.forcing.default_forcing _num_type = Val{info.helpers.numbers.num_type}() - _forcing_vars_info = info.settings.forcing.variables + _forcing_vars_info = info.experiment.data_settings.forcing.variables outcubes = mapCube(psTEMYax, (in_cubes_all...,); @@ -166,7 +166,7 @@ function run_psTEMYax(selected_models::Tuple, forcing::NamedTuple, in_cube_param clean_data=(; _data_fill, _forcing_default_info, _num_type, _forcing_vars_info), indims=indims, outdims=run_helpers.output_dims, - max_cache=info.settings.experiment.exe_rules.yax_max_cache, + max_cache=info.experiment.exe_rules.yax_max_cache, ispar=true, ) return outcubes From ce607770c5d0a93b2114fc0780f63ade7c5d0972 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Thu, 24 Jul 2025 16:58:43 +0200 Subject: [PATCH 19/27] yet another, input arg spinup_TEM --- lib/SindbadTEM/src/runTEMCube.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/SindbadTEM/src/runTEMCube.jl b/lib/SindbadTEM/src/runTEMCube.jl index 73b141581..88c95a3f4 100644 --- a/lib/SindbadTEM/src/runTEMCube.jl +++ b/lib/SindbadTEM/src/runTEMCube.jl @@ -20,8 +20,8 @@ function coreTEMYax(selected_models, loc_forcing, loc_land, tem_info) spinup_forcing = getAllSpinupForcing(loc_forcing, tem_info.spinup_sequence, tem_info); land_prec = definePrecomputeTEM(selected_models, loc_forcing_t, loc_land, tem_info.model_helpers) - land_prec = precomputeTEM(selected_models, loc_forcing_t, land_prec, tem_info.model_helpers) - land_spin = spinupTEM(selected_models, spinup_forcing, loc_forcing_t, land_prec, tem_info) # ? not good? see git diff + land_prec = precomputeTEM(selected_models, loc_forcing_t, land_prec, tem_info.model_helpers) # ? do I need this step here? + land_spin = spinupTEM(selected_models, spinup_forcing, loc_forcing_t, land_prec, tem_info, run_helpers.tem_info.run.spinup_TEM) # land_spin = land_prec land_time_series = timeLoopTEM(selected_models, loc_forcing, loc_forcing_t, land_spin, tem_info, tem_info.run.debug_model) From 8cd71fdcc357b390c17de883f399fc8001ccce4c Mon Sep 17 00:00:00 2001 From: Lazaro Date: Thu, 24 Jul 2025 17:58:57 +0200 Subject: [PATCH 20/27] set plantForm, and constants explicitly because, why not? --- examples/exp_global_runs/fire_models.jl | 3 ++- .../settings_global_runs/model_structure.json | 6 ++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/examples/exp_global_runs/fire_models.jl b/examples/exp_global_runs/fire_models.jl index 413f4e625..ba881e751 100644 --- a/examples/exp_global_runs/fire_models.jl +++ b/examples/exp_global_runs/fire_models.jl @@ -1,4 +1,5 @@ -fire_models = (:constants, +fire_models = ( + :constants, :wCycleBase, :rainSnow, :rainIntensity, diff --git a/examples/exp_global_runs/settings_global_runs/model_structure.json b/examples/exp_global_runs/settings_global_runs/model_structure.json index b903b17f4..9162be600 100644 --- a/examples/exp_global_runs/settings_global_runs/model_structure.json +++ b/examples/exp_global_runs/settings_global_runs/model_structure.json @@ -153,6 +153,9 @@ "PET": { "approach": "Lu2005" }, + "plantForm": { + "approach": "PFT" + }, "rainSnow": { "approach": "Tair" }, @@ -217,6 +220,9 @@ "wCycleBase": { "approach": "simple" }, + "constants": { + "approach": "numbers" + }, "WUE": { "approach": "expVPDDayCo2" } From 833746eedd2b83459f578154bc5bcfe4c1f65c30 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Thu, 24 Jul 2025 18:14:23 +0200 Subject: [PATCH 21/27] no run_helpers at this level --- lib/SindbadTEM/src/runTEMCube.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SindbadTEM/src/runTEMCube.jl b/lib/SindbadTEM/src/runTEMCube.jl index 88c95a3f4..31a0a6917 100644 --- a/lib/SindbadTEM/src/runTEMCube.jl +++ b/lib/SindbadTEM/src/runTEMCube.jl @@ -21,7 +21,7 @@ function coreTEMYax(selected_models, loc_forcing, loc_land, tem_info) land_prec = definePrecomputeTEM(selected_models, loc_forcing_t, loc_land, tem_info.model_helpers) land_prec = precomputeTEM(selected_models, loc_forcing_t, land_prec, tem_info.model_helpers) # ? do I need this step here? - land_spin = spinupTEM(selected_models, spinup_forcing, loc_forcing_t, land_prec, tem_info, run_helpers.tem_info.run.spinup_TEM) + land_spin = spinupTEM(selected_models, spinup_forcing, loc_forcing_t, land_prec, tem_info, tem_info.run.spinup_TEM) # land_spin = land_prec land_time_series = timeLoopTEM(selected_models, loc_forcing, loc_forcing_t, land_spin, tem_info, tem_info.run.debug_model) From 3103e1b7ffa37e4d03928eb977f784c746aad246 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Thu, 24 Jul 2025 18:49:24 +0200 Subject: [PATCH 22/27] also include WROASTED with Mortality scheme --- .../cCycleDisturbance_WROASTEDMortality.jl | 152 ++++++++++++++++++ 1 file changed, 152 insertions(+) create mode 100644 src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl diff --git a/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl b/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl new file mode 100644 index 000000000..db77cc717 --- /dev/null +++ b/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl @@ -0,0 +1,152 @@ +export cCycleDisturbance_WROASTEDMortality + +#! format: off +struct cCycleDisturbance_WROASTEDMortality <: cCycleDisturbance end +#! format: on + +function define(params::cCycleDisturbance_WROASTEDMortality, forcing, land, helpers) + @unpack_nt begin + (c_giver, c_taker) ⇐ land.constants + (cVeg, cEco) ⇐ land.pools + zix ⇐ helpers.pools + (z_zero, o_one) ⇐ land.constants + end + zix_veg_all = Tuple(vcat(getZix(cVeg, helpers.pools.zix.cVeg)...)) + c_lose_to_zix_vec = Tuple{Int}[] + for zixVeg ∈ zix_veg_all + # make reserve pool flow to slow litter pool/woody debris + if helpers.pools.components.cEco[zixVeg] == :cVegReserve + c_lose_to_zix = helpers.pools.zix.cLitSlow + else + c_lose_to_zix = c_taker[[(c_giver .== zixVeg)...]] + end + ndxNoVeg = Int[] + for ndxl ∈ c_lose_to_zix + if ndxl ∉ zix_veg_all + push!(ndxNoVeg, ndxl) + end + end + push!(c_lose_to_zix_vec, Tuple(ndxNoVeg)) + end + c_lose_to_zix_vec = Tuple(c_lose_to_zix_vec) + + # initialize disturbance outputs + # c_Veg_Mortality = zero.(cVeg) + c_Veg_Mortality = zero.(cEco) + c_Fire_Flux = zero.(cEco) + cFireTotal = z_zero + + @pack_nt begin + (zix_veg_all, c_lose_to_zix_vec) ⇒ land.cCycleDisturbance + (c_Veg_Mortality, c_Fire_Flux) ⇒ land.diagnostics + cFireTotal ⇒ land.fluxes + end + return land +end + +function compute(params::cCycleDisturbance_WROASTEDMortality, forcing, land, helpers) + ## unpack disturbance variables + @unpack_nt begin + # for vegetation die-off + c_fVegDieOff ⇐ land.diagnostics + # for fires + (c_Veg_Mortality, c_Fire_Flux, c_fire_fba) ⇐ land.diagnostics + # for fires + (c_Fire_cci, c_Fire_k) ⇐ land.diagnostics + cEco ⇐ land.pools + cFireTotal ⇐ land.fluxes + zix ⇐ helpers.pools + c_remain ⇐ land.states + (zix_veg_all, c_lose_to_zix_vec) ⇐ land.cCycleDisturbance # TODO: double check the new flow for fire, are indices correct? + (c_giver, c_taker) ⇐ land.constants + (z_zero, o_one) ⇐ land.constants + c_model ⇐ land.models + end + # set c_Fire_Flux and c_Veg_Mortality and cFireTotal to 0 + cFireTotal = z_zero + for izix in zix.cEco + @rep_elem z_zero ⇒ (c_Fire_Flux, izix, :cEco) + @rep_elem z_zero ⇒ (c_Veg_Mortality, izix, :cEco) + end + + # if there is not fire and no dieoff, pack and return + if c_fire_fba != 0.0f0 || c_fVegDieOff != 0.0f0 + # @show c_fire_fba, c_fVegDieOff + # compute vegetation mortality, and splits to litter and atmosphere + for zixVeg ∈ zix_veg_all + # total mortality fraction of vegetation pool + f_loss = c_fVegDieOff + c_fire_fba * c_Fire_k[zixVeg] + cLoss = maxZero(cEco[zixVeg] - c_remain) * f_loss + # part that is combusted and that goes to the litter pools + cLossFire = cLoss * (c_fire_fba * c_Fire_k[zixVeg]) / f_loss * c_Fire_cci[zixVeg] # ? if f_loss is zero this is undefined + cLossSoil = cLoss - cLossFire + # deplet pool + @add_to_elem -cLoss ⇒ (cEco, zixVeg, :cEco) + # transfer non combusted part + c_lose_to_zix = c_lose_to_zix_vec[zixVeg] + for tZ ∈ eachindex(c_lose_to_zix) + tarZix = c_lose_to_zix[tZ] + toGain = cLossSoil / oftype(cLossSoil, length(c_lose_to_zix)) + @add_to_elem toGain ⇒ (cEco, tarZix, :cEco) + end + # feed c_Fire_Flux and c_Veg_Mortality (@rep_elem) + @rep_elem cLossFire ⇒ (c_Fire_Flux, zixVeg, :cEco) + @rep_elem cLoss ⇒ (c_Veg_Mortality, zixVeg, :cEco) + end + + # compute fire flux from litter and soils + for zixDead ∈ (zix.cLit..., zix.cSoil...) + # total combustion from pool + f_loss = c_fire_fba * c_Fire_cci[zixDead] + cLoss = maxZero(cEco[zixDead] * f_loss) + cLossFire = cLoss + # deplet pool + @add_to_elem -cLoss ⇒ (cEco, zixDead, :cEco) # ? this one is also a new addition + # Print at every time step, left and right! + # feed c_Fire_Flux + @rep_elem cLossFire ⇒ (c_Fire_Flux, zixDead, :cEco) + end + # total fire flux + cFireTotal = totalS(c_Fire_Flux) + end + + ## pack land variables + @pack_nt begin + cEco ⇒ land.pools + cFireTotal ⇒ land.fluxes + (c_Veg_Mortality, c_Fire_Flux) ⇒ land.diagnostics + end + land = adjustPackPoolComponents(land, helpers, c_model) + return land +end + +@doc """ +move all vegetation carbon pools except reserve to respective flow target when there is disturbance + +# Parameters +$(SindbadParameters) + +--- + +# compute: +Disturb the carbon cycle pools using cCycleDisturbance_WROASTEDMortality + +*Inputs* + - land.pools.cEco: carbon pool at the end of spinup + +*Outputs* + +# update + +update pools and states in cCycleDisturbance_WROASTEDMortality + + - land.pools.cEco + +--- + +# Extended help + +*Created by* + - Nuno | nunocarvalhais +""" +cCycleDisturbance_WROASTEDMortality From 4ae96834a26de2878c782a6be715181fd7443a95 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Fri, 25 Jul 2025 12:09:13 +0200 Subject: [PATCH 23/27] added missing purpose, and finally, globally running again! --- .../settings_global_runs/model_structure.json | 2 +- .../cCycleDisturbance_WROASTEDMortality.jl | 22 +++---------------- 2 files changed, 4 insertions(+), 20 deletions(-) diff --git a/examples/exp_global_runs/settings_global_runs/model_structure.json b/examples/exp_global_runs/settings_global_runs/model_structure.json index 9162be600..b8e496518 100644 --- a/examples/exp_global_runs/settings_global_runs/model_structure.json +++ b/examples/exp_global_runs/settings_global_runs/model_structure.json @@ -63,7 +63,7 @@ "use_in_spinup": false }, "cCycleDisturbance": { - "approach": "WROASTED", + "approach": "WROASTEDMortality", "use_in_spinup": false }, "cFlow": { diff --git a/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl b/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl index db77cc717..8e95f42b8 100644 --- a/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl +++ b/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl @@ -120,27 +120,11 @@ function compute(params::cCycleDisturbance_WROASTEDMortality, forcing, land, hel return land end -@doc """ -move all vegetation carbon pools except reserve to respective flow target when there is disturbance - -# Parameters -$(SindbadParameters) - ---- - -# compute: -Disturb the carbon cycle pools using cCycleDisturbance_WROASTEDMortality +purpose(::Type{cCycleDisturbance_WROASTEDMortality}) = "This is used for vegetation die-off and fire disturbance events. Moves carbon in reserve pool to slow litter pool, and all other carbon pools except reserve pool to their respective carbon flow target pools during disturbance events." -*Inputs* - - land.pools.cEco: carbon pool at the end of spinup - -*Outputs* - -# update - -update pools and states in cCycleDisturbance_WROASTEDMortality +@doc """ - - land.pools.cEco +$(getModelDocString(cCycleDisturbance_WROASTEDMortality)) --- From 7b99f9e117e46bac7171dfec7ea5036232764bec Mon Sep 17 00:00:00 2001 From: Lazaro Date: Fri, 25 Jul 2025 16:33:13 +0200 Subject: [PATCH 24/27] slurm jobs, global runs --- .../exp_global_runs/exp_global_parameters.jl | 22 +++++- .../exp_global_runs/exp_global_parameters.sh | 9 +++ .../exp_global_runs/exp_global_pft_runs.jl | 71 +++++++++++++++++++ .../exp_global_runs/exp_global_pft_runs.sh | 23 ++++++ examples/exp_global_runs/exp_global_runs.jl | 22 +++++- examples/exp_global_runs/exp_global_runs.sh | 23 ++++++ 6 files changed, 168 insertions(+), 2 deletions(-) create mode 100644 examples/exp_global_runs/exp_global_parameters.sh create mode 100644 examples/exp_global_runs/exp_global_pft_runs.jl create mode 100644 examples/exp_global_runs/exp_global_pft_runs.sh create mode 100644 examples/exp_global_runs/exp_global_runs.sh diff --git a/examples/exp_global_runs/exp_global_parameters.jl b/examples/exp_global_runs/exp_global_parameters.jl index 85d70b298..5f7b40fdc 100644 --- a/examples/exp_global_runs/exp_global_parameters.jl +++ b/examples/exp_global_runs/exp_global_parameters.jl @@ -34,4 +34,24 @@ ps_path = "/Net/Groups/BGI/work_5/scratch/lalonso/parameters_ALL_new_0d25.zarr" out_params = mapParamsAll([cube_PFTs, cube_KG, ds_cubes_in], trainedNN, lower_bound, upper_bound, ps_names, ps_path; metadata_global= metadata_global - ) \ No newline at end of file + ) + +# TODO: run only with PFTs as covariates! and neural network trained with PFTs only +# ? load trained neural network +path_model_pft = joinpath("/Net/Groups/BGI/work_5/scratch/lalonso/checkpoint_epoch_208.jld2") #! update path! +trainedNN, lower_bound, upper_bound, ps_names, metadata_global = + loadTrainedNN(path_model_pft) + +## ? load global PFT dataset +ds = open_dataset("/Net/Groups/BGI/work_5/scratch/lalonso/CovariatesGlobal_025.zarr") +cube_PFTs = readcubedata(ds["PFT_mask"]) +cube_PFTs = Float32.(cube_PFTs) +# ? match spatial grid dimensions +cube_PFTs = YAXArray((longitude(-179.875:0.25:179.875), latitude(89.875:-0.25:-89.875)), + cube_PFTs.data, cube_PFTs.properties) + + # ? compute and save new parameters +ps_path = "/Net/Groups/BGI/work_5/scratch/lalonso/parameters_PFT_new_0d25.zarr" + +out_params = mapParamsPFT(cube_PFTs, trainedNN, lower_bound, upper_bound, ps_names, ps_path; + metadata_global= metadata_global) \ No newline at end of file diff --git a/examples/exp_global_runs/exp_global_parameters.sh b/examples/exp_global_runs/exp_global_parameters.sh new file mode 100644 index 000000000..2dd71593d --- /dev/null +++ b/examples/exp_global_runs/exp_global_parameters.sh @@ -0,0 +1,9 @@ +#!/bin/bash +#SBATCH --job-name=Params +#SBATCH --ntasks=10 # Total number of tasks +#SBATCH --cpus-per-task=10 # 10 CPUs per task +#SBATCH --mem-per-cpu=3GB # 3GB per CPU +#SBATCH --time=23:50:00 # 10 minutes runtime + +# telling slurm where to write output and error +#SBATCH -o /Net/Groups/BGI/tscratch/lalonso/SindbadOutput/Params_slurm-%A_%a.out diff --git a/examples/exp_global_runs/exp_global_pft_runs.jl b/examples/exp_global_runs/exp_global_pft_runs.jl new file mode 100644 index 000000000..f98197f18 --- /dev/null +++ b/examples/exp_global_runs/exp_global_pft_runs.jl @@ -0,0 +1,71 @@ +# ENV["JULIA_NUM_PRECOMPILE_TASKS"] = "1" +using Distributed +using SlurmClusterManager # pkg> add https://github.com/lazarusA/SlurmClusterManager.jl.git#la/asynclaunch +addprocs(SlurmManager()) + +using Sindbad +using SindbadUtils +using SindbadTEM +using SindbadData +using SindbadSetup +using SindbadML.JLD2 + +@everywhere begin + import Pkg + Pkg.activate(@__DIR__) + using Sindbad + using SindbadUtils + using SindbadTEM + using SindbadData + using SindbadSetup + using SindbadML.JLD2 +end + +@everywhere println("Worker $(myid()): Number of Threads = ", Threads.nthreads()) + +include("fire_models.jl"); + +experiment_output = "PFT_025" +remote_local = "/Net/Groups/BGI/tscratch/lalonso/" +remote_root = joinpath(remote_local, experiment_output) +mkpath(remote_root) + +domain = "Global"; +replace_info_spatial = Dict( + "experiment.model_output.path" => remote_root, + "experiment.basics.domain" => domain * "_spatial", + "experiment.flags.spinup_TEM" => true, + "experiment.flags.run_optimization" => false, + "info.settings.experiment.flags.calc_cost" => false, + "model_structure.sindbad_models" => fire_models + ); + +experiment_json = "../exp_global_runs/settings_global_runs/experiment.json"; + +info = getExperimentInfo(experiment_json; replace_info=replace_info_spatial); + +path_model = joinpath("/Net/Groups/BGI/work_5/scratch/lalonso/checkpoint_epoch_208.jld2") #! update path! +model_props = JLD2.load(path_model) +tbl_params = model_props["parameter_table"] + +forcing = getForcing(info); +# ds = open_dataset("/Net/Groups/BGI/work_4/scratch/lalonso/GlobalForcingSet.zarr"); +run_helpers = prepTEM(forcing, info); + +tuple_models = getTupleFromLongTuple(info.models.forward); + +# ? newly generated parameters with PFTs only as covariates +ps_path = "/Net/Groups/BGI/work_5/scratch/lalonso/parameters_PFT_new_0d25.zarr" + +in_cube_params = Cube(ps_path) +in_cube_ps = permutedims(in_cube_params, (2,3,1)) +in_cube_ps = readcubedata(in_cube_ps) + +yax_max_cache = info.experiment.exe_rules.yax_max_cache + +outcubes = SindbadTEM.runTEMYaxParameters( + tuple_models, + forcing, + in_cube_ps, + tbl_params, + info); \ No newline at end of file diff --git a/examples/exp_global_runs/exp_global_pft_runs.sh b/examples/exp_global_runs/exp_global_pft_runs.sh new file mode 100644 index 000000000..a37750a2e --- /dev/null +++ b/examples/exp_global_runs/exp_global_pft_runs.sh @@ -0,0 +1,23 @@ +#!/bin/bash +#SBATCH --job-name=PFT +#SBATCH --ntasks=10 # Total number of tasks +#SBATCH --cpus-per-task=10 # 10 CPUs per task +#SBATCH --mem-per-cpu=3GB # 3GB per CPU +#SBATCH --time=23:50:00 # 10 minutes runtime + +# telling slurm where to write output and error +#SBATCH -o /Net/Groups/BGI/tscratch/lalonso/SindbadOutput/PFT_slurm-%A_%a.out + +# if needed load modules here +module load proxy +module load julia/1.11 + +# if needed add export variables here +export JULIA_NUM_THREADS=${SLURM_CPUS_PER_TASK} + +################ +# +# run the program +# +################ +julia --project --heap-size-hint=16G exp_global_pft_runs.jl diff --git a/examples/exp_global_runs/exp_global_runs.jl b/examples/exp_global_runs/exp_global_runs.jl index dc38241ab..feb69a097 100644 --- a/examples/exp_global_runs/exp_global_runs.jl +++ b/examples/exp_global_runs/exp_global_runs.jl @@ -1,4 +1,11 @@ # ENV["JULIA_NUM_PRECOMPILE_TASKS"] = "1" +import Pkg +Pkg.activate(@__DIR__) + +using Distributed +using SlurmClusterManager # pkg> add https://github.com/lazarusA/SlurmClusterManager.jl.git#la/asynclaunch +addprocs(SlurmManager()) + using Sindbad using SindbadUtils using SindbadTEM @@ -6,6 +13,19 @@ using SindbadData using SindbadSetup using SindbadML.JLD2 +@everywhere begin + import Pkg + Pkg.activate(@__DIR__) + using Sindbad + using SindbadUtils + using SindbadTEM + using SindbadData + using SindbadSetup + using SindbadML.JLD2 +end + +@everywhere println("Worker $(myid()): Number of Threads = ", Threads.nthreads()) + include("fire_models.jl"); experiment_output = "ALL_025" @@ -37,7 +57,7 @@ run_helpers = prepTEM(forcing, info); tuple_models = getTupleFromLongTuple(info.models.forward); -# newly generated parameters +# ? newly generated parameters ps_path = "/Net/Groups/BGI/work_5/scratch/lalonso/parameters_ALL_new_0d25.zarr" in_cube_params = Cube(ps_path) diff --git a/examples/exp_global_runs/exp_global_runs.sh b/examples/exp_global_runs/exp_global_runs.sh new file mode 100644 index 000000000..cb1b3467e --- /dev/null +++ b/examples/exp_global_runs/exp_global_runs.sh @@ -0,0 +1,23 @@ +#!/bin/bash +#SBATCH --job-name=ALL +#SBATCH --ntasks=10 # Total number of tasks +#SBATCH --cpus-per-task=10 # 10 CPUs per task +#SBATCH --mem-per-cpu=3GB # 3GB per CPU +#SBATCH --time=23:50:00 # 10 minutes runtime + +# telling slurm where to write output and error +#SBATCH -o /Net/Groups/BGI/tscratch/lalonso/SindbadOutput/ALL_slurm-%A_%a.out + +# if needed load modules here +module load proxy +module load julia/1.11 + +# if needed add export variables here +export JULIA_NUM_THREADS=${SLURM_CPUS_PER_TASK} + +################ +# +# run the program +# +################ +julia --project --heap-size-hint=16G exp_global_runs.jl From 75dd4672bf7b72a7ef2a99c44d0cbf1be743b3ab Mon Sep 17 00:00:00 2001 From: Lazaro Date: Fri, 25 Jul 2025 18:23:47 +0200 Subject: [PATCH 25/27] vegetated land masks --- .../PlantFunctionalTypes_sinusoidal.jl | 70 +++++++++++++++++++ .../PlantFunctionalTypes_sinusoidal_map.jl | 25 +++++++ .../global_masks_area_veg/Project.toml | 4 ++ .../global_masks_area_veg/README_masks.md | 53 ++++++++++++++ .../global_masks_area_veg/VegetatedLand.jl | 44 ++++++++++++ .../VegetatedLand_map.jl | 38 ++++++++++ .../VegetatedLand_resample_4326.jl | 38 ++++++++++ .../VegetatedLand_resample_4326_map.jl | 31 ++++++++ .../global_masks_area_veg/area_mask.jl | 29 ++++++++ 9 files changed, 332 insertions(+) create mode 100644 examples/exp_global_runs/global_masks_area_veg/PlantFunctionalTypes_sinusoidal.jl create mode 100644 examples/exp_global_runs/global_masks_area_veg/PlantFunctionalTypes_sinusoidal_map.jl create mode 100644 examples/exp_global_runs/global_masks_area_veg/Project.toml create mode 100644 examples/exp_global_runs/global_masks_area_veg/README_masks.md create mode 100644 examples/exp_global_runs/global_masks_area_veg/VegetatedLand.jl create mode 100644 examples/exp_global_runs/global_masks_area_veg/VegetatedLand_map.jl create mode 100644 examples/exp_global_runs/global_masks_area_veg/VegetatedLand_resample_4326.jl create mode 100644 examples/exp_global_runs/global_masks_area_veg/VegetatedLand_resample_4326_map.jl create mode 100644 examples/exp_global_runs/global_masks_area_veg/area_mask.jl diff --git a/examples/exp_global_runs/global_masks_area_veg/PlantFunctionalTypes_sinusoidal.jl b/examples/exp_global_runs/global_masks_area_veg/PlantFunctionalTypes_sinusoidal.jl new file mode 100644 index 000000000..86baafbfe --- /dev/null +++ b/examples/exp_global_runs/global_masks_area_veg/PlantFunctionalTypes_sinusoidal.jl @@ -0,0 +1,70 @@ +using Pkg +Pkg.activate(@__DIR__) +Pkg.add(url="https://github.com/EarthyScience/UnpackSinTiles.jl", rev="la/split_land_cover_types") + +using YAXArrays, Zarr, FillArrays +using DimensionalData +using Dates +using DelimitedFiles +using ProgressMeter +using UnpackSinTiles + +# ? https://lpdaac.usgs.gov/documents/1006/MCD64_User_Guide_V61.pdf +# Appendix B Coordinate conversion for the MODIS sinusoidal projection +# Navigation of the tiled MODIS products in the sinusoidal projection can be performed using the forward +# and inverse mapping transformations described here. We’ll first need to define a few constants: +# R = 6371007.181 m, the radius of the idealized sphere representing the Earth; +# T = 1111950 m, the height and width of each MODIS tile in the projection plane; +# xmin = -20015109 m, the western limit of the projection plane; +# ymax = 10007555 m, the northern limit of the projection plane; +# w = T /2400= 463.31271653 m, the actual size of a “500-m” MODIS sinusoidal grid cell. +# ! full ellipsoid +#SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs") + + +properties = Dict("LandCover" => "MODIS/MCD64A1.061", + "crs" => "+proj=sinu +lon_0=0 +type=crs", + "grid_cell_width" => "463.31271653 m", + "R" => "6371007.181 m") + +axlist = ( + Y(range(10007555, step=-463.31271653, length=43200)), + X(range(-20015109, step=463.31271653, length=86400)), +) + +skeleton_data = Fill(NaN32, 43200, 86400) + +cube = YAXArray(axlist, skeleton_data, properties) + +path_to_lc = "/Net/Groups/BGI/work_5/scratch/lalonso/PlantFunctionalTypes.zarr" + +ds = YAXArrays.Dataset(; (:Plant_Functional_Type => cube,)...) + +d_cube = savedataset(ds; path=path_to_lc, driver=:zarr, skeleton=true, overwrite=true) + +# test one, open one! +hv_tile = "h01v08" +in_date = "2009.01.01" + +root_path = "/Net/Groups/BGI/data/DataStructureMDI/DATA/Incoming/MODIS/MCD12Q1.061/orgdata/" +lc_type1 = loadTileVariable(hv_tile, in_date, root_path, "LC_Type1") + +#! convert those 0xXY +LC_TypeX_vals = Float32.(lc_type1) + +# ? open and update entries in dataset +outar = zopen(joinpath(path_to_lc, "Plant_Functional_Type"), "w") + +mtiles = modisTiles()[:] + +@showprogress for hv_tile in mtiles + lc_type1 = loadTileVariable(hv_tile, in_date, root_path, "LC_Type1") + if !isnothing(lc_type1) + LC_TypeX_vals = Float32.(lc_type1) + # get indices in array for the current tile + hv_indices = getTileInterval(hv_tile; h_bound=1:86400, v_bound=1:43200) + #! update entries! + @info hv_tile + outar[hv_indices] = LC_TypeX_vals + end +end \ No newline at end of file diff --git a/examples/exp_global_runs/global_masks_area_veg/PlantFunctionalTypes_sinusoidal_map.jl b/examples/exp_global_runs/global_masks_area_veg/PlantFunctionalTypes_sinusoidal_map.jl new file mode 100644 index 000000000..9e20e40e5 --- /dev/null +++ b/examples/exp_global_runs/global_masks_area_veg/PlantFunctionalTypes_sinusoidal_map.jl @@ -0,0 +1,25 @@ +using CairoMakie +using YAXArrays +using Zarr +path_masks = "/Net/Groups/BGI/work_5/scratch/lalonso/maps_masks" + +# ? load plant funcional types +path_to_lc = "/Net/Groups/BGI/work_5/scratch/lalonso/PlantFunctionalTypes.zarr" + +ds_open = open_dataset(path_to_lc) +cube_open = ds_open["Plant_Functional_Type"] +cube_pft = readcubedata(cube_open) #! be careful, this is ~14G in memory + +let + using CairoMakie + fig = Figure(figure_padding=(0); size = (8640, 4320)) + ax = Axis(fig[1,1]) + hidedecorations!(ax) + hidespines!(ax) + heatmap!(ax, cube_pft[1:10:end, 1:10:end], + colormap=Categorical(cgrad(:ground_cover, 16, categorical=true)), # tol_land_cover + colorrange = (1,16), highclip=:grey95, + lowclip= :grey13 + ) + save(joinpath(path_masks, "PlantFunctionalTypes.png"), fig) +end \ No newline at end of file diff --git a/examples/exp_global_runs/global_masks_area_veg/Project.toml b/examples/exp_global_runs/global_masks_area_veg/Project.toml new file mode 100644 index 000000000..29ce34199 --- /dev/null +++ b/examples/exp_global_runs/global_masks_area_veg/Project.toml @@ -0,0 +1,4 @@ +[deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c" +Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99" diff --git a/examples/exp_global_runs/global_masks_area_veg/README_masks.md b/examples/exp_global_runs/global_masks_area_veg/README_masks.md new file mode 100644 index 000000000..4a3f0bc02 --- /dev/null +++ b/examples/exp_global_runs/global_masks_area_veg/README_masks.md @@ -0,0 +1,53 @@ +# Why this README? +Well, just to state that we have raw data in high resolution for land cover types, which we use for the PFT scenario. +And when computing burn areas we need to know how much of that is actullay vegetated, hence we should know what is the fraction of that area to use in the final calculation, this is done by going up from the highest resolution! + +This folder contains all the scripts used to pre-proccess all that! + +## UnpackSinTiles.jl + +> [!CAUTION] +> For this, I created another package to reads and stacks all tiles into global maps in a lon/lat grid. +> The original products are in the MODIS sinusoidal projection. + +For some scripts you will need to add (specially when reading the original files): + +```julia +pkg> add https://github.com/EarthyScience/UnpackSinTiles.jl#la/split_land_cover_types +``` + +it will be included where appropiate, however in this folder's `Project.toml` is not included! + +> [!INFO] +> Also, some scripts for that are already available in that package, however for completeness workflow I'm copying them here! + +### IGBP Vegetation Types Classification: `LC_Type1` + +> [!IMPORTANT] +> The user guide uses an scale from 0-16, zero being water. But the actual data goes from 1 to 17, where now 17 is water. +> See also https://www.ceom.ou.edu/static/docs/IGBP.pdf + +> [!TIP] +> Non-vegetated areas: `Dominant Plant Form` +> 17, 13, 15, 16, 255 + +| Value | Land Cover Classification | Dominant Plant Form | +|-------|--------------------------------------| -------------------- | +| 1 | Evergreen Needleleaf Forests | Tree | +| 2 | Evergreen Broadleaf Forests | Tree | +| 3 | Deciduous Needleleaf Forests | Tree | +| 4 | Deciduous Broadleaf Forests | Tree | +| 5 | Mixed Forests | Tree | +| 6 | Closed Shrublands | Shrub | +| 7 | Open Shrublands | Shrub | +| 8 | Woody Savannas | Savanna | +| 9 | Savannas | Savanna | +| 10 | Grasslands | Herb | +| 11 | Permanent Wetlands | Herb | +| 12 | Croplands | Herb | +| 13 | Urban and Built-up Lands | Non-Veg | +| 14 | Cropland/Natural Vegetation Mosaics | Herb | +| 15 | Permanent Snow and Ice | Non-Veg | +| 16 | Barren or Sparsely Vegetated | Non-Veg | +| 17 | Water | Non-Veg | +| 255 | Unclassified (No Data) | Non-Veg | \ No newline at end of file diff --git a/examples/exp_global_runs/global_masks_area_veg/VegetatedLand.jl b/examples/exp_global_runs/global_masks_area_veg/VegetatedLand.jl new file mode 100644 index 000000000..ff3cc4304 --- /dev/null +++ b/examples/exp_global_runs/global_masks_area_veg/VegetatedLand.jl @@ -0,0 +1,44 @@ +using YAXArrays, Zarr, FillArrays +using DimensionalData +using Dates +using DelimitedFiles +using ProgressMeter +using UnpackSinTiles +using Rasters + +function compute_mask_not(m, to_vals::AbstractVector) + if ismissing(m) + return false + else + return m .∉ Ref(to_vals) + end +end + +# Load data in SINUSOIDAL_CRS +path_to_lc = "/Net/Groups/BGI/work_5/scratch/lalonso/PlantFunctionalTypes.zarr" +cube_open = open_dataset(path_to_lc)["Plant_Functional_Type"] + +# cube_pft = readcubedata(cube_open) + +veg_forms = Dict{String, AbstractArray}() +veg_forms["non_veg"] = [17, 13, 15, 16, 255, NaN, 1.0f32] +veg_forms["tree"] = [1,2,3,4,5] +veg_forms["shrub"] = [6,7] +veg_forms["savanna"] = [8,9] +veg_forms["herb"] = [10,11,12,14] + +vegetated_land = "/Net/Groups/BGI/work_5/scratch/lalonso/VegetatedLand.zarr" + +r = mapCube(cube_open, indims=InDims(), outdims=OutDims(; overwrite=true, path=vegetated_land)) do xout, x + _veg_bool = ismissing(x) ? NaN : x + xout .= compute_mask_not(_veg_bool, veg_forms["non_veg"]) ? 1 : NaN +end + +veg_forms["non_veg_p16"] = [17, 13, 15, 255, NaN, 1.0f32] # include barren -> 16 (by removing from this list :D) + +vegetated_land_barren = "/Net/Groups/BGI/work_5/scratch/lalonso/VegetatedLand_pBarren.zarr" + +r2 = mapCube(cube_open, indims=InDims(), outdims=OutDims(; overwrite=true, path=vegetated_land_barren)) do xout, x + _veg_bool = ismissing(x) ? NaN : x + xout .= compute_mask_not(_veg_bool, veg_forms["non_veg_p16"]) ? 1 : NaN +end \ No newline at end of file diff --git a/examples/exp_global_runs/global_masks_area_veg/VegetatedLand_map.jl b/examples/exp_global_runs/global_masks_area_veg/VegetatedLand_map.jl new file mode 100644 index 000000000..ba31f35c2 --- /dev/null +++ b/examples/exp_global_runs/global_masks_area_veg/VegetatedLand_map.jl @@ -0,0 +1,38 @@ +using CairoMakie +using YAXArrays +using Zarr +path_masks = "/Net/Groups/BGI/work_5/scratch/lalonso/maps_masks" + +vegetated_land = "/Net/Groups/BGI/work_5/scratch/lalonso/VegetatedLand.zarr" +r = open_dataset(vegetated_land)["layer"] +cube_veg_land = readcubedata(r) #! be careful, this is ~14G in memory + +let + fig = Figure(figure_padding=(0); size = (8640, 4320)) + ax = Axis(fig[1,1]) + hidedecorations!(ax) + hidespines!(ax) + heatmap!(ax, cube_veg_land[1:10:end, 1:10:end], + colormap=Categorical(cgrad([:grey45, :black], 2, categorical=true)), # tol_land_cover + # colorrange = (1,16), highclip=:black, + # lowclip= :grey13 + ) + save(joinpath(path_masks, "VegetatedLand.png"), fig) +end + +vegetated_land_barren = "/Net/Groups/BGI/work_5/scratch/lalonso/VegetatedLand_pBarren.zarr" +r2 = open_dataset(vegetated_land_barren)["layer"] +cube_veg_land_barren = readcubedata(r2) #! be careful, this is ~14G in memory + +let + fig = Figure(figure_padding=(0); size = (8640, 4320)) + ax = Axis(fig[1,1]) + hidedecorations!(ax) + hidespines!(ax) + heatmap!(ax, cube_veg_land_barren[1:10:end, 1:10:end], + colormap=Categorical(cgrad([:grey45, :black], 2, categorical=true)), # tol_land_cover + # colorrange = (1,16), highclip=:black, + # lowclip= :grey13 + ) + save(joinpath(path_masks,"VegetatedLand_plus_barren.png"), fig) +end diff --git a/examples/exp_global_runs/global_masks_area_veg/VegetatedLand_resample_4326.jl b/examples/exp_global_runs/global_masks_area_veg/VegetatedLand_resample_4326.jl new file mode 100644 index 000000000..fa04e1abf --- /dev/null +++ b/examples/exp_global_runs/global_masks_area_veg/VegetatedLand_resample_4326.jl @@ -0,0 +1,38 @@ +using YAXArrays, Zarr +using DimensionalData +using Rasters +using Rasters.Lookups +using DimensionalData.Lookups +using ProgressMeter + +function yaxRaster(yax) + x_range = lookup(yax, :X).data + y_range = lookup(yax, :Y).data + _data = replace(yax.data, NaN=>0) + return Raster(_data, (Y(y_range; sampling=Intervals(Start())), X(x_range; sampling=Intervals(Start()))), + crs=ProjString("+proj=sinu +lon_0=0 +type=crs")) +end + + +vegetated_land = "/Net/Groups/BGI/work_5/scratch/lalonso/VegetatedLand.zarr" +ds_veg = open_dataset(vegetated_land) + +sin_ras = yaxRaster(readcubedata(ds_veg["layer"])) + +resampled = resample(sin_ras; size=(720, 1440), crs=EPSG(4326), method="average") + +locus_resampled = DimensionalData.shiftlocus(Center(), resampled) +new_dims = (lat(lookup(locus_resampled, :Y)), lon(lookup(locus_resampled, :X))) +zeros_to_nan = replace(locus_resampled, 0 => NaN32) + +properties = Dict{String, Any}() +properties["VegetatedLand"] = "Percentage of vegetated area per pixel" +properties["PRODUCT"] = "MODIS/MCD64A1.061" +properties["crs"] = "EPSG:4326" + +resampled_yax = YAXArray(new_dims, zeros_to_nan.data, properties) +ds_sampled = YAXArrays.Dataset(; (:VegLand => resampled_yax, )...) + +vegetated_land_0d25 = "/Net/Groups/BGI/work_5/scratch/lalonso/VegetatedLand_0d25.zarr" + +savedataset(ds_sampled, path=vegetated_land_0d25, driver=:zarr, overwrite=true) \ No newline at end of file diff --git a/examples/exp_global_runs/global_masks_area_veg/VegetatedLand_resample_4326_map.jl b/examples/exp_global_runs/global_masks_area_veg/VegetatedLand_resample_4326_map.jl new file mode 100644 index 000000000..d521fa577 --- /dev/null +++ b/examples/exp_global_runs/global_masks_area_veg/VegetatedLand_resample_4326_map.jl @@ -0,0 +1,31 @@ +using CairoMakie +using YAXArrays +using Zarr +path_masks = "/Net/Groups/BGI/work_5/scratch/lalonso/maps_masks" +# ? load vegetated land fraction +vegetated_land_0d25 = "/Net/Groups/BGI/work_5/scratch/lalonso/VegetatedLand_0d25.zarr" +cube_land = open_dataset(vegetated_land_0d25)["VegLand"] +cube_veg_land = readcubedata(cube_land) + +let + using CairoMakie + fig = Figure(; size=(1440, 720)) + ax = Axis(fig[1,1]) + hidedecorations!(ax) + hidespines!(ax) + hm = heatmap!(ax, -1*(cube_veg_land .- 1 .- 1e-3); colorscale=log10, colorrange=(1e-3, 1), + colormap=:linear_worb_100_25_c53_n256,) + Colorbar(fig[1,2], hm) + save(joinpath(path_masks, "vegetated_land_fraction_flipped_scale.png"), fig) +end + +let + using CairoMakie + fig = Figure(; size=(1440, 720)) + ax = Axis(fig[1,1]) + hidedecorations!(ax) + hidespines!(ax) + hm = heatmap!(ax, cube_veg_land; colormap=Reverse(:gist_stern), colorrange=(0,1)) + Colorbar(fig[1,2], hm) + save(joinpath(path_masks,"vegetated_land_fraction_0_to_1.png"), fig) +end diff --git a/examples/exp_global_runs/global_masks_area_veg/area_mask.jl b/examples/exp_global_runs/global_masks_area_veg/area_mask.jl new file mode 100644 index 000000000..b65c370b2 --- /dev/null +++ b/examples/exp_global_runs/global_masks_area_veg/area_mask.jl @@ -0,0 +1,29 @@ +# ? area_mask +using Rasters +using Rasters.Lookups +using Proj +using YAXArrays +using Zarr + +xdim = X(Projected(-180:0.25:180-0.25; sampling=Intervals(Start()), crs=EPSG(4326))) +ydim = Y(Projected(90-0.25:-0.25:-90; sampling=Intervals(Start()), crs=EPSG(4326))) +myraster = rand(xdim, ydim) +cs = cellarea(myraster) +area_mask = cs.data + +locus_area = DimensionalData.shiftlocus(Center(), cs) + +new_dims = (lon(lookup(locus_area, :X)), lat(lookup(locus_area, :Y))) + +properties = Dict{String, Any}() +properties["area"] = "Pixel area" +properties["units"] = "m2" +properties["crs"] = "EPSG:4326" + +area_yax = YAXArray(new_dims, locus_area.data, properties) + +ds_area = YAXArrays.Dataset(; (:area_mask => area_yax, )...) + +area_mask_0d25 = "/Net/Groups/BGI/work_5/scratch/lalonso/AreaMask_0d25.zarr" + +savedataset(ds_area, path=area_mask_0d25, driver=:zarr, overwrite=true) \ No newline at end of file From 8a23ff583c74a1eac8f03d3c2700d475f390f700 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Fri, 25 Jul 2025 18:56:26 +0200 Subject: [PATCH 26/27] fire events counts workflow, from high resolution to coarse --- .../Fire_events_counts/README_fire.md | 8 +++ .../Fire_events_counts/fill_skeleton_zarr.jl | 39 +++++++++++++++ .../Fire_events_counts/fill_skeleton_zarr.sh | 24 +++++++++ .../fire_reproject_to_4326.jl | 35 +++++++++++++ .../Fire_events_counts/fire_skeleton_zarr.jl | 49 +++++++++++++++++++ 5 files changed, 155 insertions(+) create mode 100644 examples/exp_global_runs/Fire_events_counts/README_fire.md create mode 100644 examples/exp_global_runs/Fire_events_counts/fill_skeleton_zarr.jl create mode 100644 examples/exp_global_runs/Fire_events_counts/fill_skeleton_zarr.sh create mode 100644 examples/exp_global_runs/Fire_events_counts/fire_reproject_to_4326.jl create mode 100644 examples/exp_global_runs/Fire_events_counts/fire_skeleton_zarr.jl diff --git a/examples/exp_global_runs/Fire_events_counts/README_fire.md b/examples/exp_global_runs/Fire_events_counts/README_fire.md new file mode 100644 index 000000000..83969838f --- /dev/null +++ b/examples/exp_global_runs/Fire_events_counts/README_fire.md @@ -0,0 +1,8 @@ +# Counting MODIS fire events + +The raw data comes in `hdf` files which we can read and aggregate with `UnpackSinTiles.jl`. The workflow breaks down into: + +- Creating a final fire skeleton zarr file of a given resolution, into which we will save the counts. See `fire_skeleton_zarr.jl`, where to `skeletons` are created. +- Load a single tile across all available dates, count, aggregate and update the zarr file. See `fill_skeleton_zarr.jl` +- Launch a slurm job for each tile, such that we do all of them in parallel. See `fill_skeleton_zarr.sh`. +- Reproject data from Sinusoidal to EPSG(4326). See `fire_reproject_to_4326.jl` \ No newline at end of file diff --git a/examples/exp_global_runs/Fire_events_counts/fill_skeleton_zarr.jl b/examples/exp_global_runs/Fire_events_counts/fill_skeleton_zarr.jl new file mode 100644 index 000000000..635635ea2 --- /dev/null +++ b/examples/exp_global_runs/Fire_events_counts/fill_skeleton_zarr.jl @@ -0,0 +1,39 @@ +using Pkg +Pkg.activate(@__DIR__) +Pkg.add(url="https://github.com/EarthyScience/UnpackSinTiles.jl", rev="la/split_land_cover_types") + +using UnpackSinTiles +using YAXArrays, Zarr +using DimensionalData +using SparseArrays +using DelimitedFiles + +indx = parse(Int, ARGS[1]) +# indx = 1 +hvs = readdlm(joinpath(@__DIR__, "all_tiles.txt"))[:,1] +hv_tile = hvs[indx] + +root_path = "/Net/Groups/BGI/data/DataStructureMDI/DATA/Incoming/MODIS/MCD64A1.061/MCD64A1.061/" + +# load all events for current tile +fire_events_bool = burnTimeSpan(2001, 2023, hv_tile, root_path; variable = "Burn Date") + +# weighted FireEvents +afterBurn = [spzeros(Float32, 2400, 2400) for i in 1:size(fire_events_bool, 1)]; +updateAfterBurn!(fire_events_bool, afterBurn; burnWindow = 30) +fire_events_bool = nothing + +# aggegates to 0.25 resolution +tsteps = size(afterBurn, 1) +μBurn = fill(NaN32, 40, 40, tsteps); +agg_μBurn = aggTile(μBurn, afterBurn, tsteps; res = 60) + +# open .zarr with `w` permissions +path_to_events = "/Net/Groups/BGI/tscratch/lalonso/MODIS_FIRE_EVENTS/MODIS_MCD64A1_FIRE_EVENTS_2001_2023.zarr" +outar = zopen(joinpath(path_to_events, "fire_frac"), "w") + +# get indices in array for the current tile +hv_indices = getTileInterval(hv_tile) + +# update entries! +outar[hv_indices, :] = agg_μBurn \ No newline at end of file diff --git a/examples/exp_global_runs/Fire_events_counts/fill_skeleton_zarr.sh b/examples/exp_global_runs/Fire_events_counts/fill_skeleton_zarr.sh new file mode 100644 index 000000000..1bf81b7b6 --- /dev/null +++ b/examples/exp_global_runs/Fire_events_counts/fill_skeleton_zarr.sh @@ -0,0 +1,24 @@ +#!/bin/bash +#SBATCH --job-name MODIS_FIRES +#SBATCH --cpus-per-task=2 +#SBATCH --mem-per-cpu=27G +#SBATCH --time=00:55:00 +#SBATCH --array=1-268 + +# telling slurm where to write output and error +#SBATCH -o /Net/Groups/BGI/tscratch/lalonso/MODIS_FIRE_EVENTS/events-%A_%a.out + +# if needed load modules here +module load proxy +module load julia + +# if needed add export variables here +export JULIA_NUM_THREADS=${SLURM_CPUS_PER_TASK} + +################ +# +# run the program +# +################ +sleep $SLURM_ARRAY_TASK_ID +julia --project --heap-size-hint=50G fill_skeleton_zarr.jl $SLURM_ARRAY_TASK_ID \ No newline at end of file diff --git a/examples/exp_global_runs/Fire_events_counts/fire_reproject_to_4326.jl b/examples/exp_global_runs/Fire_events_counts/fire_reproject_to_4326.jl new file mode 100644 index 000000000..00dc52168 --- /dev/null +++ b/examples/exp_global_runs/Fire_events_counts/fire_reproject_to_4326.jl @@ -0,0 +1,35 @@ +using YAXArrays, Zarr +using DimensionalData +using Rasters +using DimensionalData.Lookups +using ProgressMeter + +function yaxRaster(yax) + # x_range = lookup(yax, :X).data + # y_range = lookup(yax, :Y).data + x_range = LinRange(-20015109.355798, 20015109.355798, 1441)[1:end-1] + y_range = LinRange(10007554.677899, -10007554.677899, 721)[2:end] + + crs_sinu = yax.properties["crs"] + ras_out = replace(x -> x==1.0f32 ? NaN : x, yax) # "_FillValue= 1.0f32" + + new_raster = Raster(ras_out.data, (Y(y_range; sampling=Intervals(Start())), X(x_range; sampling=Intervals(Start()))), + crs=ProjString(crs_sinu)) + + return new_raster +end +# open .zarr +path_to_events = "/Net/Groups/BGI/tscratch/lalonso/MODIS_FIRE_EVENTS/MODIS_MCD64A1_FIRE_EVENTS_2001_2023.zarr" +ds = open_dataset(path_to_events) + +# save new values +path_to_latlon = "/Net/Groups/BGI/tscratch/lalonso/MODIS_FIRE_EVENTS/MODIS_MCD64A1_FIRE_EVENTS_2001_2023_LATLON.zarr" +outar = zopen(joinpath(path_to_latlon, "fire_frac"), "w") + +# calculate new values +@showprogress for t_index in 1:size(ds["fire_frac"], 3) + ras = yaxRaster(ds["fire_frac"][Ti=t_index]) + resampled = resample(ras; size=(720, 1440), crs=EPSG(4326), method="average") + + outar[:, :, t_index] = resampled.data +end \ No newline at end of file diff --git a/examples/exp_global_runs/Fire_events_counts/fire_skeleton_zarr.jl b/examples/exp_global_runs/Fire_events_counts/fire_skeleton_zarr.jl new file mode 100644 index 000000000..89d0b0d8d --- /dev/null +++ b/examples/exp_global_runs/Fire_events_counts/fire_skeleton_zarr.jl @@ -0,0 +1,49 @@ +using YAXArrays, Zarr, FillArrays +using DimensionalData +using Dates + +#SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs") +axlist = ( + Y(range(9.979756529777847e6, -1.0007111969122082e7, 720)), + X(range(-2.0015109355797417e7, 1.998725401355172e7, 1440)), + Ti(Date("2001-01-01"):Day(1):Date("2023-12-31")), +) + +skeleton_data = Fill(NaN32, 720, 1440, 8400) + +properties = Dict("FireEvents" => "MODIS/MCD64A1.061", + "crs" => "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs") + +c = YAXArray(axlist, skeleton_data, properties) + +path_to_events = "/Net/Groups/BGI/tscratch/lalonso/MODIS_FIRE_EVENTS/MODIS_MCD64A1_FIRE_EVENTS_2001_2023.zarr" + +ds = YAXArrays.Dataset(; (:fire_frac => c,)...) +d_cube = savedataset(ds; path=path_to_events, driver=:zarr, skeleton=true, overwrite=true) + +# ? note here: We should define the axis with the right properties. Well, it goes away, once saved. Report this as an issue. + +# axlist = ( +# Dim{:lat}(Projected(range(89.875, -89.875, 720), span = Regular(), crs = EPSG(4326), +# order = ReverseOrdered(), sampling = Intervals(Center()))), +# Dim{:lon}(Projected(range(-179.875, 179.875, 1440), span = Regular(), crs = EPSG(4326), +# order = ForwardOrdered(), sampling = Intervals(Center()))), +# Ti(Date("2001-01-01"):Day(1):Date("2023-12-31")), +# ) + +axlist = ( + Dim{:lat}(range(89.875, -89.875, 720)), + Dim{:lon}(range(-179.875, 179.875, 1440)), + Ti(Date("2001-01-01"):Day(1):Date("2023-12-31")), +) + +skeleton_data = Fill(NaN32, 720, 1440, 8400) + +properties = Dict("FireEvents" => "MODIS/MCD64A1.061",) + +c = YAXArray(axlist, skeleton_data, properties) + +path_to_events = "/Net/Groups/BGI/tscratch/lalonso/MODIS_FIRE_EVENTS/MODIS_MCD64A1_FIRE_EVENTS_2001_2023_LATLON.zarr" + +ds = YAXArrays.Dataset(; (:fire_frac => c,)...) +d_cube = savedataset(ds; path=path_to_events, driver=:zarr, skeleton=true, overwrite=true) \ No newline at end of file From 3cf44aaf64faa9ff65bab31f0761554a21a207a4 Mon Sep 17 00:00:00 2001 From: Lazaro Date: Fri, 25 Jul 2025 19:08:19 +0200 Subject: [PATCH 27/27] gfas baseline --- .../analysis_global/GFAS_co2_aggregate.jl | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 examples/exp_global_runs/analysis_global/GFAS_co2_aggregate.jl diff --git a/examples/exp_global_runs/analysis_global/GFAS_co2_aggregate.jl b/examples/exp_global_runs/analysis_global/GFAS_co2_aggregate.jl new file mode 100644 index 000000000..67ccdad54 --- /dev/null +++ b/examples/exp_global_runs/analysis_global/GFAS_co2_aggregate.jl @@ -0,0 +1,38 @@ +using YAXArrays +using Statistics +using DimensionalData +using Rasters +using Rasters.Lookups +# using NCDatasets +using NetCDF +using Proj +using JLD2 +using Dates + +# create area mask +xdim = X(Projected(0:0.1:359.9; sampling=Intervals(Start()), crs=EPSG(4326))) +ydim = Y(Projected(89.9:-0.1:-90; sampling=Intervals(Start()), crs=EPSG(4326))) +myraster = rand(xdim, ydim) +area_mask = cellarea(myraster) +area_mask = area_mask.data + +# ? load files! + +path_gfas = "/Net/Groups/data_BGC/gfas/0d10_daily/co2fire/2003/" +gfas_files = readdir(path_gfas) + +sum_gfas_co2_cat = [] +for _year in 2003:2022 + path_gfas = "/Net/Groups/data_BGC/gfas/0d10_daily/co2fire/$_year" + gfas_files = readdir(path_gfas) + for gfas_file in gfas_files + yax_one = Cube(joinpath(path_gfas, gfas_file)) + _sum_gfas = mapslices(x -> sum(skipmissing(x .* area_mask * 24 * 60 * 60)), yax_one, + dims=("longitude", "latitude")) + push!(sum_gfas_co2_cat, _sum_gfas) + @info gfas_file + end +end + +# sum_gfas_co2_cat = reduce(vcat, sum_gfas_co2_cat) +# jldsave("co2_global.jld2"; co2_global=sum_gfas_co2_cat) \ No newline at end of file