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" 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 diff --git a/examples/exp_global_runs/Project.toml b/examples/exp_global_runs/Project.toml new file mode 100644 index 000000000..59ae968f0 --- /dev/null +++ b/examples/exp_global_runs/Project.toml @@ -0,0 +1,13 @@ +[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" +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" +YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c" +Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99" 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 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..5f7b40fdc --- /dev/null +++ b/examples/exp_global_runs/exp_global_parameters.jl @@ -0,0 +1,57 @@ +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 = joinpath("/Net/Groups/BGI/work_5/scratch/lalonso/checkpoint_epoch_208.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_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 + ) + +# 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 new file mode 100644 index 000000000..feb69a097 --- /dev/null +++ b/examples/exp_global_runs/exp_global_runs.jl @@ -0,0 +1,74 @@ +# 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 +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" +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"] + +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 +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/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 diff --git a/examples/exp_global_runs/fire_models.jl b/examples/exp_global_runs/fire_models.jl new file mode 100644 index 000000000..ba881e751 --- /dev/null +++ b/examples/exp_global_runs/fire_models.jl @@ -0,0 +1,87 @@ +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/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 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..971bfc0db --- /dev/null +++ b/examples/exp_global_runs/settings_global_runs/experiment.json @@ -0,0 +1,90 @@ +{ + "basics": { + "config_files": { + "forcing": "forcing.json", + "model_structure": "model_structure.json", + "optimization": "optimization.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..b8e496518 --- /dev/null +++ b/examples/exp_global_runs/settings_global_runs/model_structure.json @@ -0,0 +1,264 @@ +{ + "default_model": { + "implicit_t_repeat": 1, + "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": "GSI_PlantForm" + }, + "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": "WROASTEDMortality", + "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" + }, + "plantForm": { + "approach": "PFT" + }, + "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" + }, + "constants": { + "approach": "numbers" + }, + "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 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 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/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..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) @@ -97,11 +114,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 +128,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 +197,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 eb4cdcafc..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 @@ -430,16 +433,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.experiment.exe_rules.land_output_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 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 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 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)) diff --git a/lib/SindbadSetup/src/setupOutput.jl b/lib/SindbadSetup/src/setupOutput.jl index b0553986f..29aee30d3 100644 --- a/lib/SindbadSetup/src/setupOutput.jl +++ b/lib/SindbadSetup/src/setupOutput.jl @@ -255,8 +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)") + 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) 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] 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/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, diff --git a/lib/SindbadTEM/src/runTEMCube.jl b/lib/SindbadTEM/src/runTEMCube.jl index c46a57fb4..31a0a6917 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) # ? do I need this step here? + 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) @@ -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 + + +""" + runTEMYaxParameters(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 + 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.experiment.data_settings.forcing.default_forcing + _num_type = Val{info.helpers.numbers.num_type}() + _forcing_vars_info = info.experiment.data_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.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 diff --git a/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl b/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl new file mode 100644 index 000000000..8e95f42b8 --- /dev/null +++ b/src/Models/cCycleDisturbance/cCycleDisturbance_WROASTEDMortality.jl @@ -0,0 +1,136 @@ +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 + +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." + +@doc """ + +$(getModelDocString(cCycleDisturbance_WROASTEDMortality)) + +--- + +# Extended help + +*Created by* + - Nuno | nunocarvalhais +""" +cCycleDisturbance_WROASTEDMortality diff --git a/src/Models/cFireBurnedArea/cFireBurnedArea.jl b/src/Models/cFireBurnedArea/cFireBurnedArea.jl new file mode 100644 index 000000000..90a465f9a --- /dev/null +++ b/src/Models/cFireBurnedArea/cFireBurnedArea.jl @@ -0,0 +1,18 @@ +export cFireBurnedArea + +abstract type cFireBurnedArea <: LandEcosystem end +purpose(::Type{cFireBurnedArea}) = "Disturbance of the carbon cycle pools due to fire." +includeApproaches(cFireBurnedArea, @__DIR__) + +@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..558484cc8 --- /dev/null +++ b/src/Models/cFireCombustionCompleteness/cFireCombustionCompleteness.jl @@ -0,0 +1,18 @@ +export cFireCombustionCompleteness + +abstract type cFireCombustionCompleteness <: LandEcosystem end +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. + +# 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..73245451f --- /dev/null +++ b/src/Models/cFireMortality/cFireMortality.jl @@ -0,0 +1,17 @@ +export cFireMortality + +abstract type cFireMortality <: LandEcosystem end +purpose(::Type{cFireMortality}) = "Disturbance of the carbon cycle pools due to fire." +includeApproaches(cFireMortality, @__DIR__) + +@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 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 +