run globally

.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

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"

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`

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

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

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

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)

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"

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)

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)

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

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);