From 99119f861595a694497d994647479e559d4804d5 Mon Sep 17 00:00:00 2001 From: ThanhNguyen428 Date: Tue, 22 Jul 2025 13:39:53 -0700 Subject: [PATCH 1/9] for FluxnetSimulationsExt module, wrote access functions to get simulation info + parameters for 4 specific sites with hardcoded information and any other site with mapped info --- Artifacts.toml | 16 +- ext/FluxnetSimulationsExt.jl | 11 +- ext/fluxnet_simulations/generic_site.jl | 163 +++++++++++++++++- .../get_fluxnet_metadata.jl | 121 +++++++++++++ src/Artifacts.jl | 4 + src/standalone/Vegetation/Canopy.jl | 1 - .../spatially_varying_parameters.jl | 156 +++++++++++++++++ test/integrated/fluxnet_sim.jl | 91 ++++++++++ 8 files changed, 553 insertions(+), 10 deletions(-) create mode 100644 ext/fluxnet_simulations/get_fluxnet_metadata.jl diff --git a/Artifacts.toml b/Artifacts.toml index e49c7aec8e..8e79b8b354 100644 --- a/Artifacts.toml +++ b/Artifacts.toml @@ -43,17 +43,18 @@ git-tree-sha1 = "136f1db3ed969614fb589c5250545a3ed1e8aaab" url = "https://caltech.box.com/shared/static/of1admpndmikoumtgk5j3yvt92v71awk.gz" ["clm_data_0.9x1.25"] -git-tree-sha1 = "283c62220fca8c4afe36a62348d3e9a159af2ee9" +git-tree-sha1 = "38488be5dd476890deb6fbd2b26e3e1cee2d449e" [["clm_data_0.9x1.25".download]] - sha256 = "88d5899a729de800017a74bd8a7278582c2c55c0d8730a529bae384425d70e4e" - url = "https://caltech.box.com/shared/static/mxs3l1c1dppjwy81assk8w8ofn9d70pu.gz" + sha256 = "cb5ead78868243879261d26cbb42d665bfd98135ebbebc052487d4d13ebf3a98" + url = "https://caltech.box.com/shared/static/6ms2qzyiekorksewmfyie6mbonlq5ctp.gz" ["clm_data_0.125x0.125"] -git-tree-sha1 = "6284ddefbc7937d9c1fb68fa731ff3f00b68e917" +git-tree-sha1 = "b8816f0af2c6d182111a156c3b8eed2ed83efea7" [["clm_data_0.125x0.125".download]] - sha256 = "8d2a61baad53346515f4a66c4342f8aafce37ca9520020c17f99cc4f4aa98b15" - url = "https://caltech.box.com/shared/static/uteaw1l6fg7wmwhb3ey1c0jq4r8vg0ts.gz" + sha256 = "993b5b97b6ec22ec484db4e523b65a231d736ead8b5e0759ad3abd00050bb58c" + url = "https://caltech.box.com/shared/static/fx9138ysguyg4g65l4tvibbk16er0nqo.gz" + [era5_land_forcing_data2008] git-tree-sha1 = "93d2e93f491e77cb8fba2a1b8b3946f38bde469e" [era5_land_forcing_data2008_lowres] @@ -185,3 +186,6 @@ git-tree-sha1 = "c35ba0e899040cb8153226ab751f69100f475d39" [[mizoguchi_soil_freezing_data.download]] sha256 = "0027cc080ba45ba33dc790b176ec2854353ce7dce4eae4bef72963b0dd944e0b" url = "https://caltech.box.com/shared/static/tn1bnqjmegyetw5kzd2ixq5pbnb05s3u.gz" + +[fluxnet2015] +git-tree-sha1 = "94d6eb8e3bc5fc361a43597ab4fb6941bdbe2850" \ No newline at end of file diff --git a/ext/FluxnetSimulationsExt.jl b/ext/FluxnetSimulationsExt.jl index 2ec5060846..d7a864f4ef 100644 --- a/ext/FluxnetSimulationsExt.jl +++ b/ext/FluxnetSimulationsExt.jl @@ -6,12 +6,15 @@ import ClimaUtilities.TimeManager: ITime, date using Thermodynamics using Dates using DelimitedFiles -using DocStringExtensions using Insolation import ClimaLand.Parameters as LP using ClimaLand using ClimaLand.Canopy using ClimaLand.PlantHydraulics +using ClimaLand.Domains +using ClimaCore + +using NCDatasets export prescribed_forcing_fluxnet, make_set_fluxnet_initial_conditions, get_comparison_data, @@ -19,7 +22,11 @@ export prescribed_forcing_fluxnet, get_data_dt, replace_hyphen, get_parameters, - get_domain_info + get_domain_info, + get_location + +# Get the metadata for the fluxnet sites. +include("fluxnet_simulations/get_fluxnet_metadata.jl") # Include site-specific configurations, as well as the default generic site. include("fluxnet_simulations/generic_site.jl") diff --git a/ext/fluxnet_simulations/generic_site.jl b/ext/fluxnet_simulations/generic_site.jl index e2d07c46fd..035e479e2d 100644 --- a/ext/fluxnet_simulations/generic_site.jl +++ b/ext/fluxnet_simulations/generic_site.jl @@ -2,7 +2,168 @@ # MODULE FUNCTIONS # ################################### -# TODO +""" + get_domain_info(FT; dz_bottom = FT(1.5), dz_top = FT(0.1), + nelements = 20, zmin = FT(-10), zmax = FT(0)) + +Gets and returns primary domain information for a generic Fluxnet site, using +default values corresponding to a 10m deep soil column with 20 layers, with +a resolution of 10 cm at the top of the domain and 1.5m at the bottom of the domain. +""" +function FluxnetSimulations.get_domain_info( + FT; + dz_bottom = FT(1.5), + dz_top = FT(0.1), + nelements = 20, + zmin = FT(-10), + zmax = FT(0), +) + + dz_tuple = (dz_bottom, dz_top) + + return (; dz_tuple, nelements, zmin, zmax) +end + +function FluxnetSimulations.get_location( + site_ID::String; + site_info = FluxnetSimulationsExt.get_site_info(site_ID), + lat = site_info[:lat], + long = site_info[:long], + time_offset = site_info[:time_offset], +) + atmos_h = site_info[:atmospheric_sensor_height][1] # take the first height as default + return (; time_offset, lat, long, atmos_h) +end + +""" + get_parameters(FT, site_ID, domain, pft::Pft; kwargs...) + +Gets parameters for a generic Fluxnet site based on PFT (plant functional type), +supplemented with additional autofilled values from US-MOz (Missouri Ozark) site. + +Data sources: + +Hydraulics parameters: + - Holtzman, Nataniel, et al. 2023, https://doi.org/10.1029/2023WR035481 +""" +function FluxnetSimulations.get_parameters( + FT, + site_ID, + domain, + pft::Pft = ClimaLand.Canopy.clm_dominant_pft(domain.space.surface), + ; + soil_params_gupta = ClimaLand.Soil.soil_vangenuchten_parameters( + domain.space.subsurface, + FT, + ), + soil_params_grids = ClimaLand.Soil.soil_composition_parameters( + domain.space.subsurface, + FT, + ), + soil_ν = soil_params_gupta[:ν], + soil_K_sat = soil_params_gupta[:K_sat], + soil_S_s = FT(1e-2), + soil_hydrology_cm = soil_params_gupta[:hydrology_cm], + θ_r = soil_params_gupta[:θ_r], + ν_ss_quartz = soil_params_grids[:ν_ss_quartz], + ν_ss_om = soil_params_grids[:ν_ss_om], + ν_ss_gravel = soil_params_grids[:ν_ss_gravel], + z_0m_soil = FT(0.01), + z_0b_soil = FT(0.01), + soil_ϵ = FT(0.98), + soil_albedo = ClimaLand.Soil.CLMTwoBandSoilAlbedo{FT}(; + ClimaLand.Soil.clm_soil_albedo_parameters(domain.space.surface)..., + ), + Ω = pft.parameters.Ω, + χl = pft.parameters.χl, + G_Function = CLMGFunction(χl), + α_PAR_leaf = pft.parameters.α_PAR_leaf, + λ_γ_PAR = FT(5e-7), + α_NIR_leaf = pft.parameters.α_NIR_leaf, + τ_PAR_leaf = pft.parameters.τ_PAR_leaf, + τ_NIR_leaf = pft.parameters.τ_NIR_leaf, + ϵ_canopy = pft.parameters.ϵ_canopy, + ac_canopy = pft.parameters.ac_canopy, + g1 = pft.parameters.g1, + Drel = FT(1.6), + g0 = ClimaCore.Fields.field2array( + ClimaLand.Canopy.clm_medlyn_g0(domain.space.surface), + )[1], + Vcmax25 = pft.parameters.Vcmax25, + pc = FT(-2.0e6), + sc = FT(5e-6), + SAI = FT(0), + f_root_to_shoot = pft.parameters.f_root_to_shoot, + K_sat_plant = pft.parameters.K_sat_plant, + ψ63 = pft.parameters.ψ63, + Weibull_param = FT(4), + a = FT(0.1 * 0.0098), + conductivity_model = PlantHydraulics.Weibull{FT}( + K_sat_plant, + ψ63, + Weibull_param, + ), + retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a), + plant_ν = pft.parameters.plant_ν, + plant_S_s = FT(1e-2 * 0.0098), + rooting_depth = pft.parameters.rooting_depth, + n_stem = Int64(0), + n_leaf = Int64(1), + h_leaf = ClimaCore.Fields.field2array( + ClimaLand.Canopy.clm_z_top(domain.space.surface), + )[1], + h_canopy = FluxnetSimulationsExt.get_canopy_height(site_ID), + h_stem = h_canopy - h_leaf, + z0_m = FT(0.13) * h_canopy, + z0_b = FT(0.1) * z0_m, +) + return (; + soil_ν, + soil_K_sat, + soil_S_s, + soil_hydrology_cm, + θ_r, + ν_ss_quartz, + ν_ss_om, + ν_ss_gravel, + z_0m_soil, + z_0b_soil, + soil_ϵ, + soil_albedo, + Ω, + χl, + G_Function, + α_PAR_leaf, + λ_γ_PAR, + τ_PAR_leaf, + α_NIR_leaf, + τ_NIR_leaf, + ϵ_canopy, + ac_canopy, + g1, + Drel, + g0, + Vcmax25, + SAI, + f_root_to_shoot, + K_sat_plant, + ψ63, + Weibull_param, + a, + conductivity_model, + retention_model, + plant_ν, + plant_S_s, + rooting_depth, + n_stem, + n_leaf, + h_leaf, + h_stem, + h_canopy, + z0_m, + z0_b, + ) +end ################################### # UTILITIES # diff --git a/ext/fluxnet_simulations/get_fluxnet_metadata.jl b/ext/fluxnet_simulations/get_fluxnet_metadata.jl new file mode 100644 index 0000000000..f4efb1efcc --- /dev/null +++ b/ext/fluxnet_simulations/get_fluxnet_metadata.jl @@ -0,0 +1,121 @@ +using DelimitedFiles +using ClimaLand + +function parse_array_field(field) + if field == "" || field == "NaN" + return Float64[] + elseif field isa Float64 + return isnan(field) ? Float64[] : [field] + else + return parse.(Float64, split(String(field), ";")) + end +end + + +""" + get_site_info(site_ID) + +This function retrieves the site information from an artifact for a given fluxnet site ID. + +Returns a Float64-valued namedtuple containing: +- `lat`: Latitude of the site (deg) +- `long`: Longitude of the site (deg) +- `time_offset`: Time offset from UTC (hours). Note that the ClimaLand convention is for + UTC = local_time + time_offset. +- `atmospheric_sensor_height` (Vector): Heights of atmospheric sensors (m) +""" +function get_site_info(site_ID; fluxnet2015_metadata_path = nothing) + if fluxnet2015_metadata_path === nothing + # Use the default path from ClimaLand.Artifacts + fluxnet2015_data_path = ClimaLand.Artifacts.fluxnet2015_data_path() + fluxnet2015_metadata_path = + joinpath(fluxnet2015_data_path, "metadata_DD_clean.csv") + end + raw = readdlm(fluxnet2015_metadata_path, ',', Any) + header = raw[1, :] + data = raw[2:end, :] + + # Find the row matching site_ID + row_idx = findfirst(row -> row[1] == site_ID, eachrow(data)) + if isnothing(row_idx) + error("Site ID $site_ID not found in metadata.") + end + + site_metadata = data[row_idx, :] + + varnames = + ["latitude", "longitude", "utc_offset", "atmospheric_sensor_heights"] + column_name_map = Dict( + varname => findfirst(header[:] .== varname) for varname in varnames + ) + + nan_if_empty(x) = + (x isa String && isempty(x)) ? NaN : + (x isa AbstractFloat && isnan(x)) ? NaN : x + + for (varname, column) in column_name_map + field = site_metadata[column] + if field == "" || (field isa AbstractFloat && isnan(field)) + @warn "Field $(varname) is missing for site $site_ID" + end + end + + return (; + lat = nan_if_empty(site_metadata[column_name_map["latitude"]]), + long = nan_if_empty(site_metadata[column_name_map["longitude"]]), + time_offset = -1 * nan_if_empty( + site_metadata[column_name_map["utc_offset"]], + ), + atmospheric_sensor_height = parse_array_field( + site_metadata[column_name_map["atmospheric_sensor_heights"]], + ), + ) +end + + +""" + get_canopy_height(site_ID) + +This function retrieves a canopy height from an artifact for a given fluxnet site ID. + +Returns a Float64-valued namedtuple containing: +- `canopy_height`: Canopy height of the site (m) +""" +function get_canopy_height(site_ID; fluxnet2015_metadata_path = nothing) + if fluxnet2015_metadata_path === nothing + # Use the default path from ClimaLand.Artifacts + fluxnet2015_data_path = ClimaLand.Artifacts.fluxnet2015_data_path() + fluxnet2015_metadata_path = + joinpath(fluxnet2015_data_path, "metadata_DD_clean.csv") + end + + raw = readdlm(fluxnet2015_metadata_path, ',', Any) + header = raw[1, :] + data = raw[2:end, :] + + site_info = get_site_info(site_ID; fluxnet2015_metadata_path) + + # Find the row matching site_ID + row_idx = findfirst(row -> row[1] == site_ID, eachrow(data)) + if isnothing(row_idx) + error("Site ID $site_ID not found in metadata.") + end + + site_metadata = data[row_idx, :] + + varname = "canopy_height" + + col_idx = findfirst(header[:] .== varname) + + nan_if_empty(x) = + (x isa String && isempty(x)) ? NaN : + (x isa AbstractFloat && isnan(x)) ? NaN : x + + canopy_height = site_metadata[col_idx] + if canopy_height == "" || + (canopy_height isa AbstractFloat && isnan(canopy_height)) + @warn "Field $(varname) is missing for site $site_ID" + end + + return nan_if_empty(canopy_height) +end diff --git a/src/Artifacts.jl b/src/Artifacts.jl index d8c34c4ef4..70a14913e4 100644 --- a/src/Artifacts.jl +++ b/src/Artifacts.jl @@ -277,6 +277,10 @@ function experiment_fluxnet_data_path(site_ID; context = nothing) return data_path end +function fluxnet2015_data_path(; context = nothing) + return @clima_artifact("fluxnet2015", context) +end + """ esm_snowmip_data_path(; context = nothing) diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl index d54700237a..0284ccb37d 100644 --- a/src/standalone/Vegetation/Canopy.jl +++ b/src/standalone/Vegetation/Canopy.jl @@ -52,7 +52,6 @@ using Dates include("./autotrophic_respiration.jl") include("./spatially_varying_parameters.jl") - ######################################################## # Convenience constructors for Canopy model components ######################################################## diff --git a/src/standalone/Vegetation/spatially_varying_parameters.jl b/src/standalone/Vegetation/spatially_varying_parameters.jl index aded9822e8..f0d4d10a4f 100644 --- a/src/standalone/Vegetation/spatially_varying_parameters.jl +++ b/src/standalone/Vegetation/spatially_varying_parameters.jl @@ -289,3 +289,159 @@ function clm_medlyn_g1( ) return g1 end + +""" + clm_medlyn_g0( + surface_space; + regridder_type = :InterpolationsRegridder, + extrapolation_bc = ( + Interpolations.Periodic(), + Interpolations.Flat(), + Interpolations.Flat(), + ), + lowres=use_lowres_clm(surface_space), + ) + +Reads spatially varying g0 for the canopy, from a NetCDF file +based on CLM data, and regrids it to the grid defined by the +`surface_space` of the Clima simulation. Returns a NamedTuple of ClimaCore +Fields. + +The values correspond to the value of the dominant PFT at each point. + +The NetCDF files are stored in ClimaArtifacts and more detail on their origin +is provided there. The keyword arguments `regridder_type` and `extrapolation_bc` +affect the regridding by (1) changing how we interpolate to ClimaCore points which +are not in the data, and (2) changing how extrapolate to points beyond the range of the +data. The keyword argument lowres is a flag that determines if the 0.9x1.25 or 0.125x0.125 +resolution CLM data artifact is used. If the lowres flag is not provided, the clm artifact +with the closest resolution to the surface_space is used. +""" +function clm_medlyn_g0( + surface_space; + regridder_type = :InterpolationsRegridder, + extrapolation_bc = ( + Interpolations.Periodic(), + Interpolations.Flat(), + Interpolations.Flat(), + ), + lowres = use_lowres_clm(surface_space), +) + context = ClimaComms.context(surface_space) + clm_artifact_path = Artifacts.clm_data_folder_path(; context, lowres) + + g0 = SpaceVaryingInput( + joinpath(clm_artifact_path, "vegetation_properties_map.nc"), + "medlynintercept", + surface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + return g0 +end + +""" + clm_medlyn_g0( + surface_space; + regridder_type = :InterpolationsRegridder, + extrapolation_bc = ( + Interpolations.Periodic(), + Interpolations.Flat(), + Interpolations.Flat(), + ), + lowres=use_lowres_clm(surface_space), + ) + +Reads spatially varying z_top (h_leaf) for the canopy, from a NetCDF file +based on CLM data, and regrids it to the grid defined by the +`surface_space` of the Clima simulation. Returns a NamedTuple of ClimaCore +Fields. + +The values correspond to the value of the dominant PFT at each point. + +The NetCDF files are stored in ClimaArtifacts and more detail on their origin +is provided there. The keyword arguments `regridder_type` and `extrapolation_bc` +affect the regridding by (1) changing how we interpolate to ClimaCore points which +are not in the data, and (2) changing how extrapolate to points beyond the range of the +data. The keyword argument lowres is a flag that determines if the 0.9x1.25 or 0.125x0.125 +resolution CLM data artifact is used. If the lowres flag is not provided, the clm artifact +with the closest resolution to the surface_space is used. +""" +function clm_z_top( + surface_space; + regridder_type = :InterpolationsRegridder, + extrapolation_bc = ( + Interpolations.Periodic(), + Interpolations.Flat(), + Interpolations.Flat(), + ), + lowres = use_lowres_clm(surface_space), +) + context = ClimaComms.context(surface_space) + clm_artifact_path = Artifacts.clm_data_folder_path(; context, lowres) + + z_top = SpaceVaryingInput( + joinpath(clm_artifact_path, "vegetation_properties_map.nc"), + "z_top", + surface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + return z_top +end + +""" + clm_dominant_pft( + surface_space; + regridder_type = :InterpolationsRegridder, + extrapolation_bc = ( + Interpolations.Periodic(), + Interpolations.Flat(), + Interpolations.Flat(), + ), + interpolation_method = Interpolations.Constant(), + lowres = use_lowres_clm(surface_space), + ) + +Reads spatially varying dominant PFT, from a NetCDF file +based on CLM data, and regrids it to the grid defined by the +`surface_space` of the Clima simulation. Returns a Pft() object. + +The values correspond to the value of the dominant PFT at each point. + +The NetCDF files are stored in ClimaArtifacts and more detail on their origin +is provided there. The keyword arguments `regridder_type` and `extrapolation_bc` +affect the regridding by (1) changing how we interpolate to ClimaCore points which +are not in the data, and (2) changing how extrapolate to points beyond the range of the +data. The keyword argument lowres is a flag that determines if the 0.9x1.25 or 0.125x0.125 +resolution CLM data artifact is used. If the lowres flag is not provided, the clm artifact +with the closest resolution to the surface_space is used. +""" +function clm_dominant_pft( + surface_space; + regridder_type = :InterpolationsRegridder, + extrapolation_bc = ( + Interpolations.Periodic(), + Interpolations.Flat(), + Interpolations.Flat(), + ), + interpolation_method = Interpolations.Constant(), + lowres = use_lowres_clm(surface_space), +) + context = ClimaComms.context(surface_space) + clm_artifact_path = Artifacts.clm_data_folder_path(; context, lowres) + + dominant_pft_idx = SpaceVaryingInput( + joinpath(clm_artifact_path, "dominant_PFT_map.nc"), + "dominant_PFT", + surface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + + dominant_pft_idx = + trunc(Int64, ClimaCore.Fields.field2array(dominant_pft_idx)[1]) + + dominant_pft = ClimaLand.default_pfts[dominant_pft_idx] + return dominant_pft +end diff --git a/test/integrated/fluxnet_sim.jl b/test/integrated/fluxnet_sim.jl index b99a0dada1..3fe685e03d 100644 --- a/test/integrated/fluxnet_sim.jl +++ b/test/integrated/fluxnet_sim.jl @@ -5,12 +5,26 @@ import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput using ClimaLand using ClimaLand.Canopy using ClimaLand.PlantHydraulics +using ClimaLand.Domains const FT = Float64 using DelimitedFiles import ClimaLand.FluxnetSimulations as FluxnetSimulations +""" + create_site_column(FT, lat, long) + +Creates and returns a Column domain for user-given latitude and longitude coordinates. +""" +function create_site_column(FT, lat, long, dz_tuple, nelements, zmin, zmax) + + zlim = (zmin, zmax) + longlat = (long, lat) + + return Domains.Column(; zlim, nelements, longlat, dz_tuple) +end + @testset "US-Ha1 domain info + parameters" begin site_ID = FluxnetSimulations.replace_hyphen("US-Ha1") @@ -385,3 +399,80 @@ end @test h_stem == FT(0) @test z0_m == FT(0.13) * h_canopy end + +@testset "generic site domain info + parameters" begin + site_ID = "AU-Emr" + + (; time_offset, lat, long, atmos_h) = + FluxnetSimulations.get_location(site_ID) + + @test lat == FT(-23.8587) + @test long == FT(148.4746) + @test time_offset == -10 + @test atmos_h == FT(6.7) + + # domain information + (; dz_tuple, nelements, zmin, zmax) = FluxnetSimulations.get_domain_info(FT) + + @test dz_tuple == (FT(1.5), FT(0.1)) + @test nelements == 20 + @test zmin == FT(-10) + @test zmax == FT(0) + + domain = create_site_column(FT, lat, long, dz_tuple, nelements, zmin, zmax) + + (; + soil_ν, + soil_K_sat, + soil_S_s, + soil_hydrology_cm, + θ_r, + ν_ss_quartz, + ν_ss_om, + ν_ss_gravel, + z_0m_soil, + z_0b_soil, + soil_ϵ, + soil_albedo, + Ω, + χl, + G_Function, + α_PAR_leaf, + λ_γ_PAR, + τ_PAR_leaf, + α_NIR_leaf, + τ_NIR_leaf, + ϵ_canopy, + ac_canopy, + g1, + Drel, + g0, + Vcmax25, + SAI, + f_root_to_shoot, + K_sat_plant, + ψ63, + Weibull_param, + a, + conductivity_model, + retention_model, + plant_ν, + plant_S_s, + rooting_depth, + n_stem, + n_leaf, + h_leaf, + h_stem, + h_canopy, + z0_m, + z0_b, + ) = FluxnetSimulations.get_parameters(FT, site_ID, domain; θ_r = FT(0.067)) + + @test θ_r == FT(0.067) + @test Ω == 0.75 + @test g0 == FT(100.0) + @test χl == FT(-0.3) + @test h_leaf == FT(0.5) + @test Vcmax25 == FT(2.4e-5) + +end From 0fbc3a466b3b8f557c7b7b8e89565ba8f16f9a9c Mon Sep 17 00:00:00 2001 From: ThanhNguyen428 Date: Tue, 19 Aug 2025 17:20:57 -0700 Subject: [PATCH 2/9] pain --- experiments/integrated/fluxnet/run_fluxnet.jl | 19 +++++----- ext/fluxnet_simulations/US-Ha1.jl | 33 +++++------------ ext/fluxnet_simulations/US-MOz.jl | 16 +-------- ext/fluxnet_simulations/US-NR1.jl | 26 +++----------- ext/fluxnet_simulations/US-Var.jl | 16 +-------- ext/fluxnet_simulations/generic_site.jl | 18 +++++----- .../get_fluxnet_metadata.jl | 4 +-- src/Artifacts.jl | 35 ++++++++++++++++--- src/standalone/Vegetation/pfts.jl | 20 +++++------ test/integrated/fluxnet_sim.jl | 21 ++++------- 10 files changed, 86 insertions(+), 122 deletions(-) diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl index 53ec147b10..dfacbf8fb3 100644 --- a/experiments/integrated/fluxnet/run_fluxnet.jl +++ b/experiments/integrated/fluxnet/run_fluxnet.jl @@ -41,6 +41,16 @@ site_ID_val = FluxnetSimulations.replace_hyphen(site_ID) (; time_offset, lat, long) = FluxnetSimulations.get_location(FT, Val(site_ID_val)) (; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID_val)) + +# Construct the ClimaLand domain to run the simulation on +land_domain = Column(; + zlim = (zmin, zmax), + nelements = nelements, + dz_tuple = dz_tuple, + longlat = (long, lat), +) +surface_domain = ClimaLand.Domains.obtain_surface_domain(land_domain) + (; soil_ν, soil_K_sat, @@ -89,15 +99,6 @@ site_ID_val = FluxnetSimulations.replace_hyphen(site_ID) z0_b, ) = FluxnetSimulations.get_parameters(FT, Val(site_ID_val)) -# Construct the ClimaLand domain to run the simulation on -land_domain = Column(; - zlim = (zmin, zmax), - nelements = nelements, - dz_tuple = dz_tuple, - longlat = (long, lat), -) -surface_domain = ClimaLand.Domains.obtain_surface_domain(land_domain) - # Set up the timestepping information for the simulation dt = Float64(450) # 7.5 minutes diff --git a/ext/fluxnet_simulations/US-Ha1.jl b/ext/fluxnet_simulations/US-Ha1.jl index fcc992f08c..e27db57df6 100644 --- a/ext/fluxnet_simulations/US-Ha1.jl +++ b/ext/fluxnet_simulations/US-Ha1.jl @@ -32,6 +32,9 @@ end Returns geographical information for US-Ha1 (Massachusetts Harvard Forest) Fluxnet site. The values are provided as defaults, and can be overwritten by passing the corresponding keyword arguments to this function. + +Atmosphere height: + - https://atmos.seas.harvard.edu/research-harvard_forest-instrumentation """ function FluxnetSimulations.get_location( FT, @@ -39,28 +42,9 @@ function FluxnetSimulations.get_location( time_offset = 5, lat = FT(42.5378), long = FT(-72.1715), -) - return (; time_offset, lat, long) -end - -""" - get_fluxtower_height(FT, ::Val{:US_Ha1}; kwargs...) - -Returns atmosphere height for US-Ha1 (Massachusetts Harvard Forest) Fluxnet site. -The values are provided as defaults, and can be overwritten by passing the -corresponding keyword arguments to this function. - -Data sources: - -Atmosphere height: - - https://atmos.seas.harvard.edu/research-harvard_forest-instrumentation -""" -function FluxnetSimulations.get_fluxtower_height( - FT, - ::Val{:US_Ha1}; atmos_h = FT(30), ) - return (; atmos_h,) + return (; time_offset, lat, long, atmos_h) end """ @@ -85,8 +69,7 @@ function FluxnetSimulations.get_parameters( soil_ν = FT(0.5), soil_K_sat = FT(4e-7), soil_S_s = FT(1e-3), - soil_vg_n = FT(2.05), - soil_vg_α = FT(0.04), + hydrology_cm = vanGenuchten{FT}(; α = FT(0.04), n = FT(2.05)), θ_r = FT(0.067), ν_ss_quartz = FT(0.1), ν_ss_om = FT(0.1), @@ -94,8 +77,10 @@ function FluxnetSimulations.get_parameters( z_0m_soil = FT(0.01), z_0b_soil = FT(0.001), soil_ϵ = FT(0.98), - soil_α_PAR = FT(0.2), - soil_α_NIR = FT(0.2), + soil_albedo = Soil.ConstantTwoBandSoilAlbedo{FT}(; + PAR_albedo = FT(0.2), + NIR_albedo = FT(0.2), + ), Ω = FT(0.69), χl = FT(0.5), G_Function = ConstantGFunction(χl), diff --git a/ext/fluxnet_simulations/US-MOz.jl b/ext/fluxnet_simulations/US-MOz.jl index a3e3f37557..94c35bda4d 100644 --- a/ext/fluxnet_simulations/US-MOz.jl +++ b/ext/fluxnet_simulations/US-MOz.jl @@ -39,23 +39,9 @@ function FluxnetSimulations.get_location( time_offset = 7, lat = FT(38.7441), long = FT(-92.2000), -) - return (; time_offset, lat, long) -end - -""" - get_fluxtower_height(FT, ::Val{:US_MOz}; kwargs...) - -Returns atmosphere height for US-Ha1 (Missouri Ozark) Fluxnet site. -The values are provided as defaults, and can be overwritten by passing the -corresponding keyword arguments to this function. -""" -function FluxnetSimulations.get_fluxtower_height( - FT, - ::Val{:US_MOz}; atmos_h = FT(32), ) - return (; atmos_h,) + return (; time_offset, lat, long, atmos_h) end """ diff --git a/ext/fluxnet_simulations/US-NR1.jl b/ext/fluxnet_simulations/US-NR1.jl index b07e168672..5fa383a048 100644 --- a/ext/fluxnet_simulations/US-NR1.jl +++ b/ext/fluxnet_simulations/US-NR1.jl @@ -32,36 +32,20 @@ end Returns geographical information for US-NR1 (Colorado Niwot Ridge) Fluxnet site. The values are provided as defaults, and can be overwritten by passing the corresponding keyword arguments to this function. -""" -function FluxnetSimulations.get_location( - FT, - ::Val{:US_NR1}; - time_offset = 7, - lat = FT(40.0329), - long = FT(-105.5464), -) - return (; time_offset, lat, long) -end - -""" - get_fluxtower_height(FT, ::Val{:US_NR1}; kwargs...) - -Returns atmosphere height for US-NR1 (Colorado Niwot Ridge) Fluxnet site. -The values are provided as defaults, and can be overwritten by passing the -corresponding keyword arguments to this function. - -Data sources: Atmosphere height: - Metzger, Stefan & Burba, George & Burns, Sean & Blanken, Peter & Li, Jiahong & Luo, Hongyan & Zulueta, Rommel. (2016). https://doi.org/10.5194/amt-9-1341-2016 """ -function FluxnetSimulations.get_fluxtower_height( +function FluxnetSimulations.get_location( FT, ::Val{:US_NR1}; + time_offset = 7, + lat = FT(40.0329), + long = FT(-105.5464), atmos_h = FT(21.5), ) - return (; atmos_h,) + return (; time_offset, lat, long, atmos_h) end """ diff --git a/ext/fluxnet_simulations/US-Var.jl b/ext/fluxnet_simulations/US-Var.jl index b802257c25..e779b17e05 100644 --- a/ext/fluxnet_simulations/US-Var.jl +++ b/ext/fluxnet_simulations/US-Var.jl @@ -38,23 +38,9 @@ function FluxnetSimulations.get_location( time_offset = 8, lat = FT(38.4133), long = FT(-120.9508), -) - return (; time_offset, lat, long) -end - -""" - get_fluxtower_height(FT, ::Val{:US_Var}; kwargs...) - -Returns atmosphere height for US-Var (California Vaira Ranch Ione) Fluxnet site. -The values are provided as defaults, and can be overwritten by passing the -corresponding keyword arguments to this function. -""" -function FluxnetSimulations.get_fluxtower_height( - FT, - ::Val{:US_Var}; atmos_h = FT(2), ) - return (; atmos_h,) + return (; time_offset, lat, long, atmos_h) end """ diff --git a/ext/fluxnet_simulations/generic_site.jl b/ext/fluxnet_simulations/generic_site.jl index 035e479e2d..caec8a363e 100644 --- a/ext/fluxnet_simulations/generic_site.jl +++ b/ext/fluxnet_simulations/generic_site.jl @@ -60,19 +60,21 @@ function FluxnetSimulations.get_parameters( domain.space.subsurface, FT, ), - soil_ν = soil_params_gupta[:ν], - soil_K_sat = soil_params_gupta[:K_sat], + soil_ν = ClimaCore.Fields.field2array(soil_params_gupta[:ν])[1], + soil_K_sat = ClimaCore.Fields.field2array(soil_params_gupta[:K_sat])[1], soil_S_s = FT(1e-2), - soil_hydrology_cm = soil_params_gupta[:hydrology_cm], - θ_r = soil_params_gupta[:θ_r], - ν_ss_quartz = soil_params_grids[:ν_ss_quartz], - ν_ss_om = soil_params_grids[:ν_ss_om], - ν_ss_gravel = soil_params_grids[:ν_ss_gravel], + soil_hydrology_cm = ClimaCore.MatrixFields.column_field2array(soil_params_gupta[:hydrology_cm])[1], + θ_r = ClimaCore.Fields.field2array(soil_params_gupta[:θ_r])[1], + ν_ss_quartz = ClimaCore.Fields.field2array(soil_params_grids[:ν_ss_quartz])[1], + ν_ss_om = ClimaCore.Fields.field2array(soil_params_grids[:ν_ss_om])[1], + ν_ss_gravel = ClimaCore.Fields.field2array(soil_params_grids[:ν_ss_gravel])[1], z_0m_soil = FT(0.01), z_0b_soil = FT(0.01), soil_ϵ = FT(0.98), + soil_albedo_params = ClimaLand.Soil.clm_soil_albedo_parameters(domain.space.surface), soil_albedo = ClimaLand.Soil.CLMTwoBandSoilAlbedo{FT}(; - ClimaLand.Soil.clm_soil_albedo_parameters(domain.space.surface)..., + NamedTuple{keys(soil_albedo_params)}( + (ClimaCore.Fields.field2array(v)[1] for v in values(soil_albedo_params)))..., ), Ω = pft.parameters.Ω, χl = pft.parameters.χl, diff --git a/ext/fluxnet_simulations/get_fluxnet_metadata.jl b/ext/fluxnet_simulations/get_fluxnet_metadata.jl index f4efb1efcc..84a2a0bb53 100644 --- a/ext/fluxnet_simulations/get_fluxnet_metadata.jl +++ b/ext/fluxnet_simulations/get_fluxnet_metadata.jl @@ -63,9 +63,9 @@ function get_site_info(site_ID; fluxnet2015_metadata_path = nothing) return (; lat = nan_if_empty(site_metadata[column_name_map["latitude"]]), long = nan_if_empty(site_metadata[column_name_map["longitude"]]), - time_offset = -1 * nan_if_empty( + time_offset = Int64(-1 * nan_if_empty( site_metadata[column_name_map["utc_offset"]], - ), + )), atmospheric_sensor_height = parse_array_field( site_metadata[column_name_map["atmospheric_sensor_heights"]], ), diff --git a/src/Artifacts.jl b/src/Artifacts.jl index 70a14913e4..5b718f9a6b 100644 --- a/src/Artifacts.jl +++ b/src/Artifacts.jl @@ -270,11 +270,38 @@ Citation: Siyan Ma, Liukang Xu, Joseph Verfaillie, Dennis Baldocchi (2023), Amer AmeriFlux CC-BY-4.0 License """ function experiment_fluxnet_data_path(site_ID; context = nothing) - @assert site_ID ∈ ("US-MOz", "US-Var", "US-NR1", "US-Ha1") + try + full_fluxnet_path = @clima_artifact("fluxnet2015", context) + dirs = filter(d -> isdir(joinpath(full_fluxnet_path, d)), readdir(full_fluxnet_path)) + + match = findfirst(d -> occursin(site_ID, d), dirs) + @assert match !== nothing "No Fluxnet data found for site ID: $site_ID" + + site_dir = dirs[match] + parts = split(site_dir, "FULLSET") + site_path_hh = join([parts[1], "FULLSET_HH", parts[2]], "") * ".csv" + site_path_hr = join([parts[1], "FULLSET_HR", parts[2]], "") * ".csv" + + if isfile(joinpath(full_fluxnet_path, site_dir, site_path_hh)) + data_path = joinpath(full_fluxnet_path, site_dir, site_path_hh) + elseif isfile(joinpath(full_fluxnet_path, site_dir, site_path_hr)) + data_path = joinpath(full_fluxnet_path, site_dir, site_path_hr) + else + error("There exists a directory $site_dir for site ID $site_ID, but found no data files for half-hourly or hourly data.") + end + + return data_path + catch + @info "Either the full fluxnet2015 dataset does not exist locally, or the site ID was not found. + Falling back to the fluxnet_sites artifact which contains only US-MOz, US-Var, US-NR1, and US-Ha1." + + @assert site_ID ∈ ("US-MOz", "US-Var", "US-NR1", "US-Ha1") - folder_path = @clima_artifact("fluxnet_sites", context) - data_path = joinpath(folder_path, "$(site_ID).csv") - return data_path + folder_path = @clima_artifact("fluxnet_sites", context) + data_path = joinpath(folder_path, "$(site_ID).csv") + return data_path + end + end function fluxnet2015_data_path(; context = nothing) diff --git a/src/standalone/Vegetation/pfts.jl b/src/standalone/Vegetation/pfts.jl index 14364067ed..568a825534 100644 --- a/src/standalone/Vegetation/pfts.jl +++ b/src/standalone/Vegetation/pfts.jl @@ -98,7 +98,7 @@ NET_Temp = Pft( τ_NIR_leaf = 0.10, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.982, # Unitless ⋅ Ma et al. 2019 χl = 0.01, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 2234, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 2234.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 74.31, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 @@ -126,7 +126,7 @@ NET_Bor = Pft( τ_NIR_leaf = 0.10, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.982, # Unitless ⋅ Ma et al. 2019 χl = 0.01, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 2792, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 2792.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 74.31, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 @@ -154,7 +154,7 @@ NDT_Bor = Pft( τ_NIR_leaf = 0.10, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.985, # Unitless ⋅ Ma et al. 2019 χl = 0.01, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 2792, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 2792.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 74.31, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 @@ -294,7 +294,7 @@ BDT_Bor = Pft( τ_NIR_leaf = 0.25, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.968, # Unitless ⋅ Ma et al. 2019 χl = 0.25, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 745, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 745.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 140.72, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 @@ -322,7 +322,7 @@ BES_Temp = Pft( τ_NIR_leaf = 0.10, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.987, # Unitless ⋅ Ma et al. 2019 χl = 0.01, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 745, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 745.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 148.63, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 @@ -350,7 +350,7 @@ BDS_Temp = Pft( τ_NIR_leaf = 0.25, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.987, # Unitless ⋅ Ma et al. 2019 χl = 0.25, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 745, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 745.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 148.63, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 @@ -378,7 +378,7 @@ BDS_Bor = Pft( τ_NIR_leaf = 0.25, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.987, # Unitless ⋅ Ma et al. 2019 χl = 0.25, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 745, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 745.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 148.63, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 @@ -406,7 +406,7 @@ C3G_A = Pft( τ_NIR_leaf = 0.34, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.978, # Unitless ⋅ Ma et al. 2019 χl = -0.3, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 745, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 745.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 70.2, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 @@ -434,7 +434,7 @@ C3G_NA = Pft( τ_NIR_leaf = 0.34, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.978, # Unitless ⋅ Ma et al. 2019 χl = -0.3, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 745, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 745.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 166.02, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 @@ -462,7 +462,7 @@ C4G = Pft( τ_NIR_leaf = 0.34, # Unitless ⋅ CLM5.0 Table 2.3.1 ϵ_canopy = 0.978, # Unitless ⋅ Ma et al. 2019 χl = -0.3, # Unitless ⋅ CLM2 PFT Data - ac_canopy = 745, # Jm^-2K^-1 ⋅ Bonan et al. 2018 + ac_canopy = 745.0, # Jm^-2K^-1 ⋅ Bonan et al. 2018 # Stomatal Conductance Model g1 = 51.22, # kPa^(1/2) ⋅ CLM5.0 Table 2.9.1 diff --git a/test/integrated/fluxnet_sim.jl b/test/integrated/fluxnet_sim.jl index 3fe685e03d..866702debd 100644 --- a/test/integrated/fluxnet_sim.jl +++ b/test/integrated/fluxnet_sim.jl @@ -38,16 +38,14 @@ end @test zmax == FT(0) # geographical info - (; time_offset, lat, long) = + (; time_offset, lat, long, atmos_h) = FluxnetSimulations.get_location(FT, Val(site_ID)) @test time_offset == 5 @test lat == FT(42.5378) @test long == FT(-72.1715) - - (; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID)) @test atmos_h == FT(30) - + # parameters (; soil_ν, @@ -160,13 +158,11 @@ end @test zmax == FT(0) # geographical info - (; time_offset, lat, long) = + (; time_offset, lat, long, atmos_h) = FluxnetSimulations.get_location(FT, Val(site_ID), time_offset = 5) @test time_offset == 5 @test lat == FT(38.7441) @test long == FT(-92.2000) - - (; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID)) @test atmos_h == FT(32) # parameters @@ -244,14 +240,12 @@ end @test zmin == FT(-10) @test zmax == FT(0) - (; time_offset, lat, long) = + (; time_offset, lat, long, atmos_h) = FluxnetSimulations.get_location(FT, Val(site_ID)) @test time_offset == 7 @test lat == FT(40.0329) @test long == FT(-105.5464) - - (; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID)) @test atmos_h == FT(21.5) # parameters @@ -328,14 +322,13 @@ end @test nelements == 14 @test zmin == FT(-0.5) @test zmax == FT(0) + - (; time_offset, lat, long) = + (; time_offset, lat, long, atmos_h) = FluxnetSimulations.get_location(FT, Val(site_ID)) @test time_offset == 8 @test lat == FT(38.4133) @test long == FT(-120.9508) - - (; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID)) @test atmos_h == FT(2) # parameters @@ -474,5 +467,5 @@ end @test χl == FT(-0.3) @test h_leaf == FT(0.5) @test Vcmax25 == FT(2.4e-5) - end +nothing \ No newline at end of file From 6db099a01af4ce89cb539713f6cbc61ec3e2ff93 Mon Sep 17 00:00:00 2001 From: ThanhNguyen428 Date: Wed, 20 Aug 2025 13:01:16 -0700 Subject: [PATCH 3/9] run from fluxnet2015 --- experiments/integrated/fluxnet/run_fluxnet.jl | 170 ++++++++++++------ ext/FluxnetSimulationsExt.jl | 1 + ext/fluxnet_simulations/US-Ha1.jl | 8 +- ext/fluxnet_simulations/US-MOz.jl | 15 +- ext/fluxnet_simulations/US-NR1.jl | 15 +- ext/fluxnet_simulations/US-Var.jl | 15 +- ext/fluxnet_simulations/generic_site.jl | 2 +- src/Artifacts.jl | 4 +- 8 files changed, 142 insertions(+), 88 deletions(-) diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl index dfacbf8fb3..3abcde3516 100644 --- a/experiments/integrated/fluxnet/run_fluxnet.jl +++ b/experiments/integrated/fluxnet/run_fluxnet.jl @@ -33,14 +33,69 @@ if length(ARGS) < 1 end site_ID = ARGS[1] -site_ID_val = FluxnetSimulations.replace_hyphen(site_ID) -# Get the default values for this site's domain, location, and parameters -(; dz_tuple, nelements, zmin, zmax) = +if (site_ID ∈ ("US-MOz", "US-Var", "US-NR1", "US-Ha1")) + site_ID_val = FluxnetSimulations.replace_hyphen(site_ID) + + # Get the default values for this site's domain, location, and parameters + (; dz_tuple, nelements, zmin, zmax) = FluxnetSimulations.get_domain_info(FT, Val(site_ID_val)) -(; time_offset, lat, long) = + (; time_offset, lat, long, atmos_h) = FluxnetSimulations.get_location(FT, Val(site_ID_val)) -(; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID_val)) + + (; + soil_ν, + soil_K_sat, + soil_S_s, + soil_hydrology_cm, + θ_r, + ν_ss_quartz, + ν_ss_om, + ν_ss_gravel, + z_0m_soil, + z_0b_soil, + soil_ϵ, + soil_albedo, + Ω, + χl, + G_Function, + α_PAR_leaf, + λ_γ_PAR, + τ_PAR_leaf, + α_NIR_leaf, + τ_NIR_leaf, + ϵ_canopy, + ac_canopy, + g1, + Drel, + g0, + Vcmax25, + SAI, + f_root_to_shoot, + K_sat_plant, + ψ63, + Weibull_param, + a, + conductivity_model, + retention_model, + plant_ν, + plant_S_s, + rooting_depth, + n_stem, + n_leaf, + h_leaf, + h_stem, + h_canopy, + z0_m, + z0_b, + ) = FluxnetSimulations.get_parameters(FT, Val(site_ID_val)) +else + (; dz_tuple, nelements, zmin, zmax) = + FluxnetSimulations.get_domain_info(FT) + (; time_offset, lat, long, atmos_h) = + FluxnetSimulations.get_location(site_ID) +end + # Construct the ClimaLand domain to run the simulation on land_domain = Column(; @@ -51,53 +106,54 @@ land_domain = Column(; ) surface_domain = ClimaLand.Domains.obtain_surface_domain(land_domain) -(; - soil_ν, - soil_K_sat, - soil_S_s, - soil_vg_n, - soil_vg_α, - θ_r, - ν_ss_quartz, - ν_ss_om, - ν_ss_gravel, - z_0m_soil, - z_0b_soil, - soil_ϵ, - soil_α_PAR, - soil_α_NIR, - Ω, - χl, - α_PAR_leaf, - λ_γ_PAR, - τ_PAR_leaf, - α_NIR_leaf, - τ_NIR_leaf, - ϵ_canopy, - ac_canopy, - g1, - Drel, - g0, - Vcmax25, - SAI, - f_root_to_shoot, - K_sat_plant, - ψ63, - Weibull_param, - a, - conductivity_model, - retention_model, - plant_ν, - plant_S_s, - rooting_depth, - n_stem, - n_leaf, - h_leaf, - h_stem, - h_canopy, - z0_m, - z0_b, -) = FluxnetSimulations.get_parameters(FT, Val(site_ID_val)) +if (site_ID ∉ ("US-MOz", "US-Var", "US-NR1", "US-Ha1")) + (; + soil_ν, + soil_K_sat, + soil_S_s, + soil_hydrology_cm, + θ_r, + ν_ss_quartz, + ν_ss_om, + ν_ss_gravel, + z_0m_soil, + z_0b_soil, + soil_ϵ, + soil_albedo, + Ω, + χl, + G_Function, + α_PAR_leaf, + λ_γ_PAR, + τ_PAR_leaf, + α_NIR_leaf, + τ_NIR_leaf, + ϵ_canopy, + ac_canopy, + g1, + Drel, + g0, + Vcmax25, + SAI, + f_root_to_shoot, + K_sat_plant, + ψ63, + Weibull_param, + a, + conductivity_model, + retention_model, + plant_ν, + plant_S_s, + rooting_depth, + n_stem, + n_leaf, + h_leaf, + h_stem, + h_canopy, + z0_m, + z0_b, + ) = FluxnetSimulations.get_parameters(FT, site_ID, land_domain) +end # Set up the timestepping information for the simulation dt = Float64(450) # 7.5 minutes @@ -121,17 +177,17 @@ dt = Float64(450) # 7.5 minutes # Now we set up the model. For the soil model, we pick # a model type and package up parameters. soil_domain = land_domain -soil_albedo = Soil.ConstantTwoBandSoilAlbedo{FT}(; - PAR_albedo = soil_α_PAR, - NIR_albedo = soil_α_NIR, -) +# soil_albedo = Soil.ConstantTwoBandSoilAlbedo{FT}(; +# PAR_albedo = soil_α_PAR, +# NIR_albedo = soil_α_NIR, +# ) forcing = (; atmos, radiation) retention_parameters = (; ν = soil_ν, θ_r, K_sat = soil_K_sat, - hydrology_cm = vanGenuchten{FT}(; α = soil_vg_α, n = soil_vg_n), + hydrology_cm = soil_hydrology_cm, ) composition_parameters = (; ν_ss_om, ν_ss_quartz, ν_ss_gravel) diff --git a/ext/FluxnetSimulationsExt.jl b/ext/FluxnetSimulationsExt.jl index d7a864f4ef..58f7223a00 100644 --- a/ext/FluxnetSimulationsExt.jl +++ b/ext/FluxnetSimulationsExt.jl @@ -9,6 +9,7 @@ using DelimitedFiles using Insolation import ClimaLand.Parameters as LP using ClimaLand +using ClimaLand.Soil using ClimaLand.Canopy using ClimaLand.PlantHydraulics using ClimaLand.Domains diff --git a/ext/fluxnet_simulations/US-Ha1.jl b/ext/fluxnet_simulations/US-Ha1.jl index e27db57df6..837b079bb5 100644 --- a/ext/fluxnet_simulations/US-Ha1.jl +++ b/ext/fluxnet_simulations/US-Ha1.jl @@ -69,7 +69,7 @@ function FluxnetSimulations.get_parameters( soil_ν = FT(0.5), soil_K_sat = FT(4e-7), soil_S_s = FT(1e-3), - hydrology_cm = vanGenuchten{FT}(; α = FT(0.04), n = FT(2.05)), + soil_hydrology_cm = vanGenuchten{FT}(; α = FT(0.04), n = FT(2.05)), θ_r = FT(0.067), ν_ss_quartz = FT(0.1), ν_ss_om = FT(0.1), @@ -122,8 +122,7 @@ function FluxnetSimulations.get_parameters( soil_ν, soil_K_sat, soil_S_s, - soil_vg_n, - soil_vg_α, + soil_hydrology_cm, θ_r, ν_ss_quartz, ν_ss_om, @@ -131,8 +130,7 @@ function FluxnetSimulations.get_parameters( z_0m_soil, z_0b_soil, soil_ϵ, - soil_α_PAR, - soil_α_NIR, + soil_albedo, Ω, χl, G_Function, diff --git a/ext/fluxnet_simulations/US-MOz.jl b/ext/fluxnet_simulations/US-MOz.jl index 94c35bda4d..c8b67079b6 100644 --- a/ext/fluxnet_simulations/US-MOz.jl +++ b/ext/fluxnet_simulations/US-MOz.jl @@ -64,8 +64,7 @@ function FluxnetSimulations.get_parameters( soil_ν = FT(0.55), soil_K_sat = FT(4e-7), soil_S_s = FT(1e-2), - soil_vg_n = FT(2.0), - soil_vg_α = FT(0.05), + soil_hydrology_cm = vanGenuchten{FT}(; α = FT(2.0), n = FT(0.05)), θ_r = FT(0.04), ν_ss_quartz = FT(0.1), ν_ss_om = FT(0.1), @@ -73,8 +72,10 @@ function FluxnetSimulations.get_parameters( z_0m_soil = FT(0.01), z_0b_soil = FT(0.01), soil_ϵ = FT(0.98), - soil_α_PAR = FT(0.2), - soil_α_NIR = FT(0.2), + soil_albedo = Soil.ConstantTwoBandSoilAlbedo{FT}(; + PAR_albedo = FT(0.2), + NIR_albedo = FT(0.2), + ), Ω = FT(0.69), χl = FT(0.1), G_Function = CLMGFunction(χl), @@ -118,8 +119,7 @@ function FluxnetSimulations.get_parameters( soil_ν, soil_K_sat, soil_S_s, - soil_vg_n, - soil_vg_α, + soil_hydrology_cm, θ_r, ν_ss_quartz, ν_ss_om, @@ -127,8 +127,7 @@ function FluxnetSimulations.get_parameters( z_0m_soil, z_0b_soil, soil_ϵ, - soil_α_PAR, - soil_α_NIR, + soil_albedo, Ω, χl, G_Function, diff --git a/ext/fluxnet_simulations/US-NR1.jl b/ext/fluxnet_simulations/US-NR1.jl index 5fa383a048..556c74ec5d 100644 --- a/ext/fluxnet_simulations/US-NR1.jl +++ b/ext/fluxnet_simulations/US-NR1.jl @@ -71,8 +71,7 @@ function FluxnetSimulations.get_parameters( soil_ν = FT(0.45), soil_K_sat = FT(4e-7), soil_S_s = FT(1e-3), - soil_vg_n = FT(2.05), - soil_vg_α = FT(0.04), + soil_hydrology_cm = vanGenuchten{FT}(; α = FT(2.05), n = FT(0.04)), θ_r = FT(0.0), ν_ss_quartz = FT(0.1), ν_ss_om = FT(0.1), @@ -80,8 +79,10 @@ function FluxnetSimulations.get_parameters( z_0m_soil = FT(0.1), z_0b_soil = FT(0.1), soil_ϵ = FT(0.98), - soil_α_PAR = FT(0.2), - soil_α_NIR = FT(0.2), + soil_albedo = Soil.ConstantTwoBandSoilAlbedo{FT}(; + PAR_albedo = FT(0.2), + NIR_albedo = FT(0.2), + ), Ω = FT(0.71), χl = FT(0.5), G_Function = ConstantGFunction(χl), @@ -123,8 +124,7 @@ function FluxnetSimulations.get_parameters( soil_ν, soil_K_sat, soil_S_s, - soil_vg_n, - soil_vg_α, + soil_hydrology_cm, θ_r, ν_ss_quartz, ν_ss_om, @@ -132,8 +132,7 @@ function FluxnetSimulations.get_parameters( z_0m_soil, z_0b_soil, soil_ϵ, - soil_α_PAR, - soil_α_NIR, + soil_albedo, Ω, χl, G_Function, diff --git a/ext/fluxnet_simulations/US-Var.jl b/ext/fluxnet_simulations/US-Var.jl index e779b17e05..5fbc76c4f6 100644 --- a/ext/fluxnet_simulations/US-Var.jl +++ b/ext/fluxnet_simulations/US-Var.jl @@ -70,8 +70,7 @@ function FluxnetSimulations.get_parameters( soil_ν = FT(0.45), soil_K_sat = FT(0.45 / 3600 / 100), soil_S_s = FT(1e-3), - soil_vg_n = FT(2.0), - soil_vg_α = FT(2.0), + soil_hydrology_cm = vanGenuchten{FT}(; α = FT(2.0), n = FT(2.0)), θ_r = FT(0.067), ν_ss_quartz = FT(0.3), ν_ss_om = FT(0.02), @@ -79,8 +78,10 @@ function FluxnetSimulations.get_parameters( z_0m_soil = FT(0.01), z_0b_soil = FT(0.01), soil_ϵ = FT(0.98), - soil_α_PAR = FT(0.2), - soil_α_NIR = FT(0.2), + soil_albedo = Soil.ConstantTwoBandSoilAlbedo{FT}(; + PAR_albedo = FT(0.2), + NIR_albedo = FT(0.2), + ), Ω = FT(1.0), χl = FT(0.5), G_Function = ConstantGFunction(χl), @@ -124,8 +125,7 @@ function FluxnetSimulations.get_parameters( soil_ν, soil_K_sat, soil_S_s, - soil_vg_n, - soil_vg_α, + soil_hydrology_cm, θ_r, ν_ss_quartz, ν_ss_om, @@ -133,8 +133,7 @@ function FluxnetSimulations.get_parameters( z_0m_soil, z_0b_soil, soil_ϵ, - soil_α_PAR, - soil_α_NIR, + soil_albedo, Ω, χl, G_Function, diff --git a/ext/fluxnet_simulations/generic_site.jl b/ext/fluxnet_simulations/generic_site.jl index caec8a363e..fa61bfab98 100644 --- a/ext/fluxnet_simulations/generic_site.jl +++ b/ext/fluxnet_simulations/generic_site.jl @@ -115,7 +115,7 @@ function FluxnetSimulations.get_parameters( ClimaLand.Canopy.clm_z_top(domain.space.surface), )[1], h_canopy = FluxnetSimulationsExt.get_canopy_height(site_ID), - h_stem = h_canopy - h_leaf, + h_stem = ((h_canopy - h_leaf)) > 0 ? h_canopy - h_leaf : FT(0.0), z0_m = FT(0.13) * h_canopy, z0_b = FT(0.1) * z0_m, ) diff --git a/src/Artifacts.jl b/src/Artifacts.jl index 5b718f9a6b..6109661f23 100644 --- a/src/Artifacts.jl +++ b/src/Artifacts.jl @@ -232,7 +232,9 @@ end ) Return the path to the file that contains a year of fluxnet data -corresponding to a `site_ID` in the set (US-MOz, US-NR1, US-Ha1, US-Var). +corresponding to a `site_ID` in the set (US-MOz, US-NR1, US-Ha1, US-Var). If +the site_ID is not in this set, the function will try to find the corresponding +dataset, or throw an error if not found. Site publications and licenses: US-Ha1: From 4c89569e5fa9755b9578eead976d1650f30bad65 Mon Sep 17 00:00:00 2001 From: ThanhNguyen428 Date: Wed, 20 Aug 2025 13:15:10 -0700 Subject: [PATCH 4/9] testing --- .../tutorials/calibration/test_calibration.jl | 404 ++++++++++++++++++ 1 file changed, 404 insertions(+) create mode 100644 docs/src/tutorials/calibration/test_calibration.jl diff --git a/docs/src/tutorials/calibration/test_calibration.jl b/docs/src/tutorials/calibration/test_calibration.jl new file mode 100644 index 0000000000..acee6dc9ed --- /dev/null +++ b/docs/src/tutorials/calibration/test_calibration.jl @@ -0,0 +1,404 @@ +# # Model Site-Level Calibration Tutorial Using Observations +# +# In this tutorial we will calibrate the Vcmax25 and g1 parameters using latent heat +# flux observations from the FLUXNET site (US-MOz). +# +# ## Overview +# +# The tutorial covers: +# 1. Setting up a land surface model for a FLUXNET site (US-MOz) +# 2. Obtaining the observation dataset +# 3. Implementing Ensemble Kalman Inversion to calibrate Vcmax25 and g1 +# 4. Analyzing the calibration results +# +#nb # ## Prerequisites +#nb # +#nb # First, ensure you have the required packages installed: + +#nb ENV["JULIA_PKG_PRECOMPILE_AUTO"]=0 +#nb using Pkg +#nb required_pkgs = ["ClimaLand", "ClimaDiagnostics", "CairoMakie", +#nb "EnsembleKalmanProcesses", "Random", "Logging"] +#nb Pkg.add(required_pkgs) +#nb plotting_pkgs = ["ClimaAnalysis", "GeoMakie", "Printf", "StatsBase"] +#nb Pkg.add(plotting_pkgs) +#nb Pkg.instantiate() +#nb Pkg.precompile() + +# ## Setup and Imports +# +# Load all the necessary packages for land surface modeling, diagnostics, +# plotting, and ensemble methods: + +using ClimaLand +using ClimaLand.Domains: Column +using ClimaLand.Canopy +using ClimaLand.Simulations +import ClimaLand.FluxnetSimulations as FluxnetSimulations +import ClimaLand.Parameters as LP +import ClimaLand.LandSimVis as LandSimVis +import ClimaDiagnostics +import ClimaUtilities.TimeManager: date +import EnsembleKalmanProcesses as EKP +import EnsembleKalmanProcesses.ParameterDistributions as PD +using CairoMakie +CairoMakie.activate!() +using Statistics +using Logging +import Random +using Dates +using ClimaAnalysis, GeoMakie, Printf, StatsBase + +# ## Configuration and Site Setup +# +# Configure the experiment parameters and set up the FLUXNET site (US-MOz) with +# its specific location, time settings, and atmospheric conditions. + +# Set random seed for reproducibility and floating point precision +rng_seed = 1234 +rng = Random.MersenneTwister(rng_seed) +const FT = Float32 + +# Initialize land parameters and site configuration. +earth_param_set = LP.LandParameters(FT) +site_ID = "US-MOz" +site_ID_val = FluxnetSimulations.replace_hyphen(site_ID); + +# Get site-specific information: location coordinates, time offset, and sensor +# height. +(; time_offset, lat, long) = + FluxnetSimulations.get_location(FT, Val(site_ID_val)) +(; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID_val)); + +# Get maximum simulation start and end dates in UTC; these must be included in the forcing data range +start_date = DateTime(2010, 3, 1) +stop_date = DateTime(2010, 7, 1) +Δt = 450.0; # seconds + +# ## Domain and Forcing Setup +# +# Create the computational domain and load the necessary forcing data for the +# land surface model. +# ClimaLand includes the domain, forcing, and LAI as +# part of the model, but here we will need to recreate a model +# many times (for each parameter value tried), while the domain, forcing, +# and LAI are held fixed. Because of that, we define them here, once, outside of +# the function that is called to make the model. + +# Create a column domain representing a 2-meter deep soil column with 10 +# vertical layers. +zmin = FT(-2) # 2m depth +zmax = FT(0) # surface +domain = Column(; zlim = (zmin, zmax), nelements = 10, longlat = (long, lat)); + +# Load prescribed atmospheric and radiative forcing from FLUXNET data +forcing = FluxnetSimulations.prescribed_forcing_fluxnet( + site_ID, + lat, + long, + time_offset, + atmos_h, + start_date, + earth_param_set, + FT, +); + +# Get Leaf Area Index (LAI) data from MODIS satellite observations. +LAI = + ClimaLand.prescribed_lai_modis(domain.space.surface, start_date, stop_date); + +# ## Model Setup +# +# Create an integrated land model that couples canopy, snow, soil, and soil CO2 +# components. This comprehensive model allows us to simulate the full land +# surface system and its interactions. +function model(Vcmax25, g1) + Vcmax25 = FT(Vcmax25) + g1 = FT(g1) + + #md # Set up models; note: we are not using the default soil, which relies on global + #md # maps of parameters to estimate the parameters at the site. + #md # Instead we use parameter more tailored to this site. + prognostic_land_components = (:canopy, :snow, :soil, :soilco2) + #md # Create soil model + retention_parameters = (; + ν = FT(0.55), + θ_r = FT(0.04), + K_sat = FT(4e-7), + hydrology_cm = ClimaLand.Soil.vanGenuchten{FT}(; + α = FT(0.05), + n = FT(2.0), + ), + ) + composition_parameters = + (; ν_ss_om = FT(0.1), ν_ss_quartz = FT(0.1), ν_ss_gravel = FT(0.0)) + soil = ClimaLand.Soil.EnergyHydrology{FT}( + domain, + forcing, + earth_param_set; + prognostic_land_components, + additional_sources = (ClimaLand.RootExtraction{FT}(),), + runoff = ClimaLand.Soil.Runoff.SurfaceRunoff(), + retention_parameters, + composition_parameters, + S_s = FT(1e-2), + ) + surface_domain = ClimaLand.Domains.obtain_surface_domain(domain) + canopy_forcing = (; + atmos = forcing.atmos, + radiation = forcing.radiation, + ground = ClimaLand.PrognosticGroundConditions{FT}(), + ) + #md # Set up photosynthesis using the Farquhar model + photosyn_defaults = + Canopy.clm_photosynthesis_parameters(surface_domain.space.surface) + photosynthesis = Canopy.FarquharModel{FT}( + surface_domain; + photosynthesis_parameters = (; + is_c3 = photosyn_defaults.is_c3, + Vcmax25, + ), + ) + #md # Set up stomatal conductance using the Medlyn model + conductance = Canopy.MedlynConductanceModel{FT}(surface_domain; g1) + + #md # Create canopy model + canopy = Canopy.CanopyModel{FT}( + surface_domain, + canopy_forcing, + LAI, + earth_param_set; + photosynthesis, + conductance, + prognostic_land_components, + ) + + #md # Create integrated land model + land_model = + LandModel{FT}(forcing, LAI, earth_param_set, domain, Δt; soil, canopy) + + #md # Set initial conditions from FLUXNET data + set_ic! = FluxnetSimulations.make_set_fluxnet_initial_conditions( + site_ID, + start_date, + time_offset, + land_model, + ) + + #md # Configure diagnostics to output sensible and latent heat fluxes hourly + output_vars = ["lhf"] + diagnostics = ClimaLand.default_diagnostics( + land_model, + start_date; + output_writer = ClimaDiagnostics.Writers.DictWriter(), + output_vars, + average_period = :hourly, + ) + + #md # Create and run the simulation + data_dt = Second(FluxnetSimulations.get_data_dt(site_ID)) + updateat = Array(start_date:Second(Δt):stop_date) + #md # Create and run the simulation + simulation = Simulations.LandSimulation( + start_date, + stop_date, + Δt, + land_model; + set_ic!, + updateat, + user_callbacks = (), + diagnostics, + ) + solve!(simulation) + return simulation +end; + +# ## Observation and Helper Functions +# +# Define the observation function `G` that maps from parameter space to +# observation space, along with supporting functions for data processing: + +# This function runs the model and computes diurnal average of latent heat flux +function G(Vcmax25, g1) + simulation = model(Vcmax25, g1) + lhf = get_lhf(simulation) + observation = + Float64.( + get_diurnal_average( + lhf, + simulation.start_date, + simulation.start_date + Day(20), + stop_date, + ) + ) + return observation +end; + +# Helper function: Extract latent heat flux from simulation diagnostics +function get_lhf(simulation) + return ClimaLand.Diagnostics.diagnostic_as_vectors( + simulation.diagnostics[1].output_writer, + "lhf_1h_average", + ) +end; + +# Helper function: Compute diurnal average of a variable +function get_diurnal_average(var, start_date, spinup_date, stop_date) + (times, data) = var + model_dates = if times isa Vector{DateTime} + times + else + date.(times) + end + spinup_idx = findfirst(spinup_date .<= model_dates) + stop_idx = findlast(model_dates .< stop_date) + model_dates = model_dates[spinup_idx:stop_idx] + data = data[spinup_idx:stop_idx] + + hour_of_day = Hour.(model_dates) + mean_by_hour = [mean(data[hour_of_day .== Hour(i)]) for i in 0:23] + return mean_by_hour +end; + +# ## Experiment Setup +# +# We obtain observations from the FLUXNET site. The dataset contains multiple +# variables, but we will just use latent heat flux in this calibration. +dataset = FluxnetSimulations.get_comparison_data(site_ID, time_offset) +observations = get_diurnal_average( + (dataset.UTC_datetime, dataset.lhf), + start_date, + start_date + Day(20), + stop_date, +); + +# Define observation noise covariance for the ensemble Kalman process. A flat +# covariance matrix is used here for simplicity. +noise_covariance = 0.05 * EKP.I; + +# ## Prior Distribution and Calibration Configuration +# +# Set up the prior distribution for the parameter and configure the ensemble +# Kalman inversion: + +# Constrained Gaussian prior for Vcmax25 with bounds [0, 1e-3], and g1 with +# bounds [0, 1000]. The first two arguments of each prior are the mean and standard +# deviation of the Gaussian function. + +priors = [ + PD.constrained_gaussian("Vcmax25", 1e-4, 5e-5, 0, 1e-3), + PD.constrained_gaussian("g1", 150, 90, 0, 1000), +] +prior = PD.combine_distributions(priors); + +# Set the ensemble size and number of iterations +ensemble_size = 10 +N_iterations = 4; + +# ## Ensemble Kalman Inversion +# +# Initialize and run the ensemble Kalman process: + +# Sample the initial parameter ensemble from the prior distribution +initial_ensemble = EKP.construct_initial_ensemble(rng, prior, ensemble_size) + +ensemble_kalman_process = EKP.EnsembleKalmanProcess( + initial_ensemble, + observations, + noise_covariance, + EKP.Inversion(); + scheduler = EKP.DataMisfitController( + terminate_at = Inf, + on_terminate = "continue", + ), + rng, +); + +# Run the ensemble of forward models to iteratively update the parameter +# ensemble +Logging.with_logger(SimpleLogger(devnull, Logging.Error)) do + for i in 1:N_iterations + println("Iteration $i") + params_i = EKP.get_ϕ_final(prior, ensemble_kalman_process) + G_ens = hcat([G(params_i[:, j]...) for j in 1:ensemble_size]...) + EKP.update_ensemble!(ensemble_kalman_process, G_ens) + end +end; + +# ## Results Analysis and Visualization +# Get the mean of the final parameter ensemble: +EKP.get_ϕ_mean_final(prior, ensemble_kalman_process); + +# Now, let's analyze the calibration results by examining parameter evolution +# and comparing model outputs across iterations. +# +# Plot the parameter ensemble evolution over iterations to visualize +# convergence: + +dim_size = sum(length.(EKP.batch(prior))) +fig = CairoMakie.Figure(size = ((dim_size + 1) * 500, 500)) + +for i in 1:dim_size + EKP.Visualize.plot_ϕ_over_iters( + fig[1, i], + ensemble_kalman_process, + prior, + i, + ) +end + +EKP.Visualize.plot_error_over_iters( + fig[1, dim_size + 1], + ensemble_kalman_process, +) +CairoMakie.save("constrained_params_and_error.png", fig) +fig + +# Compare the model output between the first and last iterations to assess +# improvement: +fig = CairoMakie.Figure(size = (900, 400)) + +first_G_ensemble = EKP.get_g(ensemble_kalman_process, 1) +last_iter = EKP.get_N_iterations(ensemble_kalman_process) +last_G_ensemble = EKP.get_g(ensemble_kalman_process, last_iter) +n_ens = EKP.get_N_ens(ensemble_kalman_process) + +ax = Axis( + fig[1, 1]; + title = "G ensemble: first vs last iteration (n = $(n_ens), iters 1 vs $(last_iter))", + xlabel = "Observation index", + ylabel = "G", +) + +# Plot model output of first vs last iteration ensemble +for g in eachcol(first_G_ensemble) + lines!(ax, 1:length(g), g; color = (:red, 0.6), linewidth = 1.5) +end + +for g in eachcol(last_G_ensemble) + lines!(ax, 1:length(g), g; color = (:blue, 0.6), linewidth = 1.5) +end + +lines!( + ax, + 1:length(observations), + observations; + color = (:black, 0.6), + linewidth = 3, +) + +axislegend( + ax, + [ + LineElement(color = :red, linewidth = 2), + LineElement(color = :blue, linewidth = 2), + LineElement(color = :black, linewidth = 4), + ], + ["First ensemble", "Last ensemble", "Observations"]; + position = :rb, + framevisible = false, +) + +CairoMakie.resize_to_layout!(fig) +CairoMakie.save("G_first_and_last.png", fig) +fig \ No newline at end of file From 21b3e979711d8f87f6be63aee1a96aff8226f4d4 Mon Sep 17 00:00:00 2001 From: ThanhNguyen428 Date: Thu, 21 Aug 2025 11:00:44 -0700 Subject: [PATCH 5/9] calibration testing --- .../tutorials/calibration/test_calibration.jl | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/docs/src/tutorials/calibration/test_calibration.jl b/docs/src/tutorials/calibration/test_calibration.jl index acee6dc9ed..f93af2ced4 100644 --- a/docs/src/tutorials/calibration/test_calibration.jl +++ b/docs/src/tutorials/calibration/test_calibration.jl @@ -57,22 +57,21 @@ using ClimaAnalysis, GeoMakie, Printf, StatsBase # Set random seed for reproducibility and floating point precision rng_seed = 1234 rng = Random.MersenneTwister(rng_seed) -const FT = Float32 +const FT = Float64 # Initialize land parameters and site configuration. earth_param_set = LP.LandParameters(FT) -site_ID = "US-MOz" -site_ID_val = FluxnetSimulations.replace_hyphen(site_ID); +site_ID = "AU-Rob" +# site_ID_val = FluxnetSimulations.replace_hyphen(site_ID); # Get site-specific information: location coordinates, time offset, and sensor # height. -(; time_offset, lat, long) = - FluxnetSimulations.get_location(FT, Val(site_ID_val)) -(; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID_val)); +(; time_offset, lat, long, atmos_h) = + FluxnetSimulations.get_location(site_ID) # Get maximum simulation start and end dates in UTC; these must be included in the forcing data range -start_date = DateTime(2010, 3, 1) -stop_date = DateTime(2010, 7, 1) +start_date = DateTime(2014, 3, 1) +stop_date = DateTime(2014, 7, 1) Δt = 450.0; # seconds # ## Domain and Forcing Setup @@ -351,7 +350,7 @@ EKP.Visualize.plot_error_over_iters( fig[1, dim_size + 1], ensemble_kalman_process, ) -CairoMakie.save("constrained_params_and_error.png", fig) +CairoMakie.save("experiments/integrated/fluxnet/$(site_ID)/constrained_params_and_error.png", fig) fig # Compare the model output between the first and last iterations to assess @@ -400,5 +399,5 @@ axislegend( ) CairoMakie.resize_to_layout!(fig) -CairoMakie.save("G_first_and_last.png", fig) +CairoMakie.save("experiments/integrated/fluxnet/$(site_ID)/G_first_and_last.png", fig) fig \ No newline at end of file From 7e92d7db771c3a479b868a6640d224a1735e349f Mon Sep 17 00:00:00 2001 From: ThanhNguyen428 Date: Thu, 21 Aug 2025 17:08:45 -0700 Subject: [PATCH 6/9] moved files to experiments/ --- .../calibration/calibrate_all_sites.jl | 31 +++++++++++++++++ .../calibration/calibrate_all_sites.jl | 33 +++++++++++++++++++ .../calibration/test_calibration.jl | 17 ++++++---- 3 files changed, 74 insertions(+), 7 deletions(-) create mode 100644 docs/src/tutorials/calibration/calibrate_all_sites.jl create mode 100644 experiments/calibration/calibrate_all_sites.jl rename {docs/src/tutorials => experiments}/calibration/test_calibration.jl (97%) diff --git a/docs/src/tutorials/calibration/calibrate_all_sites.jl b/docs/src/tutorials/calibration/calibrate_all_sites.jl new file mode 100644 index 0000000000..30492fa6c7 --- /dev/null +++ b/docs/src/tutorials/calibration/calibrate_all_sites.jl @@ -0,0 +1,31 @@ +using ClimaLand +using ClimaLand.Artifacts +using DelimitedFiles + +# folder = "path/to/folder" +script = "test_calibration.jl" + +fluxnet2015_data_path = ClimaLand.Artifacts.fluxnet2015_data_path() +fluxnet2015_metadata_path = + joinpath(fluxnet2015_data_path, "metadata_DD_clean.csv") + + +data = readdlm(fluxnet2015_metadata_path, ',', header=true) + +header = data[2] # column names +table = data[1] # the actual rows + +site_id_col = findfirst(==("SITE_ID"), header) +site_ids = table[:, site_id_col] # vector of all site IDs +unique_site_ids = unique(site_ids) + +for f in readdir(folder; join=true) + try + println("Running $script with file: $f") + run(`julia $script $f`) + println("✅ Success on $f") + catch e + println("❌ Failed on $f: $e") + continue + end +end \ No newline at end of file diff --git a/experiments/calibration/calibrate_all_sites.jl b/experiments/calibration/calibrate_all_sites.jl new file mode 100644 index 0000000000..c8525bae53 --- /dev/null +++ b/experiments/calibration/calibrate_all_sites.jl @@ -0,0 +1,33 @@ +using ClimaLand +using ClimaLand.Artifacts +using DelimitedFiles + +script = joinpath(pkgdir(ClimaLand), "experiments/calibration/test_calibration.jl") + +# Get all site IDs from the metadata file +fluxnet2015_data_path = ClimaLand.Artifacts.fluxnet2015_data_path() +fluxnet2015_metadata_path = + joinpath(fluxnet2015_data_path, "metadata_DD_clean.csv") + +data = readdlm(fluxnet2015_metadata_path, ',', header=true) + +header = data[2] # column names +table = data[1] # the actual rows + +site_id_col = findfirst(==("site_id"), header) +site_ids = table[:, site_id_col] # vector of all site IDs +unique_site_ids = unique(site_ids) + + + +# Run the script for each unique site ID +for site_ID in unique_site_ids + println("Running $script with file: $site_ID") + try + run(`julia --project=.buildkite $script $site_ID`) + println("Success on $site_ID") + catch e + println("Failed on $site_ID: $e") + continue + end +end \ No newline at end of file diff --git a/docs/src/tutorials/calibration/test_calibration.jl b/experiments/calibration/test_calibration.jl similarity index 97% rename from docs/src/tutorials/calibration/test_calibration.jl rename to experiments/calibration/test_calibration.jl index f93af2ced4..8e66283346 100644 --- a/docs/src/tutorials/calibration/test_calibration.jl +++ b/experiments/calibration/test_calibration.jl @@ -61,7 +61,7 @@ const FT = Float64 # Initialize land parameters and site configuration. earth_param_set = LP.LandParameters(FT) -site_ID = "AU-Rob" +site_ID = ARGS[1] # site_ID_val = FluxnetSimulations.replace_hyphen(site_ID); # Get site-specific information: location coordinates, time offset, and sensor @@ -70,8 +70,8 @@ site_ID = "AU-Rob" FluxnetSimulations.get_location(site_ID) # Get maximum simulation start and end dates in UTC; these must be included in the forcing data range -start_date = DateTime(2014, 3, 1) -stop_date = DateTime(2014, 7, 1) +start_date = DateTime(2011, 3, 1) +stop_date = DateTime(2011, 7, 1) Δt = 450.0; # seconds # ## Domain and Forcing Setup @@ -331,6 +331,9 @@ EKP.get_ϕ_mean_final(prior, ensemble_kalman_process); # Now, let's analyze the calibration results by examining parameter evolution # and comparing model outputs across iterations. # +savedir = + joinpath(pkgdir(ClimaLand), "experiments/calibration/out/$(site_ID)/out") +mkpath(savedir) # Plot the parameter ensemble evolution over iterations to visualize # convergence: @@ -350,8 +353,8 @@ EKP.Visualize.plot_error_over_iters( fig[1, dim_size + 1], ensemble_kalman_process, ) -CairoMakie.save("experiments/integrated/fluxnet/$(site_ID)/constrained_params_and_error.png", fig) -fig +CairoMakie.save("$(savedir)/constrained_params_and_error.png", fig) +# fig # Compare the model output between the first and last iterations to assess # improvement: @@ -399,5 +402,5 @@ axislegend( ) CairoMakie.resize_to_layout!(fig) -CairoMakie.save("experiments/integrated/fluxnet/$(site_ID)/G_first_and_last.png", fig) -fig \ No newline at end of file +CairoMakie.save("$(savedir)/G_first_and_last.png", fig) +# fig \ No newline at end of file From 1ed1ed95e627617205ee1dc48cba4c5e1dabacdc Mon Sep 17 00:00:00 2001 From: ThanhNguyen428 Date: Fri, 5 Sep 2025 15:16:06 -0700 Subject: [PATCH 7/9] creating output files for calibration results --- .../calibration/calibrate_all_sites.jl | 1 + experiments/calibration/test_calibration.jl | 70 ++++++++++++++++--- 2 files changed, 62 insertions(+), 9 deletions(-) diff --git a/experiments/calibration/calibrate_all_sites.jl b/experiments/calibration/calibrate_all_sites.jl index c8525bae53..5612dd5c76 100644 --- a/experiments/calibration/calibrate_all_sites.jl +++ b/experiments/calibration/calibrate_all_sites.jl @@ -1,6 +1,7 @@ using ClimaLand using ClimaLand.Artifacts using DelimitedFiles +using Logging script = joinpath(pkgdir(ClimaLand), "experiments/calibration/test_calibration.jl") diff --git a/experiments/calibration/test_calibration.jl b/experiments/calibration/test_calibration.jl index 8e66283346..6b776efbeb 100644 --- a/experiments/calibration/test_calibration.jl +++ b/experiments/calibration/test_calibration.jl @@ -49,6 +49,9 @@ import Random using Dates using ClimaAnalysis, GeoMakie, Printf, StatsBase +using CSV +using DataFrames +using DelimitedFiles # ## Configuration and Site Setup # # Configure the experiment parameters and set up the FLUXNET site (US-MOz) with @@ -88,7 +91,7 @@ stop_date = DateTime(2011, 7, 1) # vertical layers. zmin = FT(-2) # 2m depth zmax = FT(0) # surface -domain = Column(; zlim = (zmin, zmax), nelements = 10, longlat = (long, lat)); +domain = Column(; zlim = (zmin, zmax), nelements = 20, longlat = (long, lat)); # Load prescribed atmospheric and radiative forcing from FLUXNET data forcing = FluxnetSimulations.prescribed_forcing_fluxnet( @@ -326,7 +329,11 @@ end; # ## Results Analysis and Visualization # Get the mean of the final parameter ensemble: -EKP.get_ϕ_mean_final(prior, ensemble_kalman_process); +final_mean_params = EKP.get_ϕ_mean_final(prior, ensemble_kalman_process); + +ϕ = EKP.get_ϕ_final(prior, ensemble_kalman_process); +y = EKP.get_obs(ensemble_kalman_process); +metrics = EKP.get_error_metrics(ensemble_kalman_process); # Now, let's analyze the calibration results by examining parameter evolution # and comparing model outputs across iterations. @@ -334,6 +341,30 @@ EKP.get_ϕ_mean_final(prior, ensemble_kalman_process); savedir = joinpath(pkgdir(ClimaLand), "experiments/calibration/out/$(site_ID)/out") mkpath(savedir) + +# Save the final parameter ensemble and error metrics to CSV files for further analysis +# or record-keeping. +writedlm( + joinpath(savedir, "final_mean_parameters_$(site_ID).csv"), + final_mean_params, + ',', +) + +writedlm( + joinpath(savedir, "final_parameter_ensemble_$(site_ID).csv"), + ϕ, + ',', +) + +writedlm( + joinpath(savedir, "latest_data_vector_$(site_ID).csv"), + y, + ',', +) + +df = DataFrame(metrics) +CSV.write(savedir, "error_metrics_$(site_ID).csv", df) + # Plot the parameter ensemble evolution over iterations to visualize # convergence: @@ -349,6 +380,7 @@ for i in 1:dim_size ) end +# TODO save numerical error over iterations EKP.Visualize.plot_error_over_iters( fig[1, dim_size + 1], ensemble_kalman_process, @@ -365,6 +397,18 @@ last_iter = EKP.get_N_iterations(ensemble_kalman_process) last_G_ensemble = EKP.get_g(ensemble_kalman_process, last_iter) n_ens = EKP.get_N_ens(ensemble_kalman_process) +writedlm( + joinpath(savedir, "first_iteration_G_ensemble_$(site_ID).csv"), + first_G_ensemble, + ',', +) + +writedlm( + joinpath(savedir, "last_iteration_G_ensemble_$(site_ID).csv"), + last_G_ensemble, + ',', +) + ax = Axis( fig[1, 1]; title = "G ensemble: first vs last iteration (n = $(n_ens), iters 1 vs $(last_iter))", @@ -381,13 +425,21 @@ for g in eachcol(last_G_ensemble) lines!(ax, 1:length(g), g; color = (:blue, 0.6), linewidth = 1.5) end -lines!( - ax, - 1:length(observations), - observations; - color = (:black, 0.6), - linewidth = 3, -) +# --- Add shaded noise around observations --- +xvals = 1:length(observations) + +# Say you have obs variance as `obs_var` (same length as `observations`) +# If you don’t, but only have noise covariance = 0.05*I, you can just use sqrt(0.05) +obs_std = sqrt.(0.05) .* length(observations) # or your own per-point std + +band!(ax, xvals, + observations .- 2 .* obs_std, + observations .+ 2 .* obs_std; + color = (:gray, 0.3), + transparency = true) + +# now lines on top +lines!(ax, xvals, observations; color = (:black, 0.8), linewidth = 3) axislegend( ax, From 4e0d28608a4da77bb2713b16824a6b3e792f2dfa Mon Sep 17 00:00:00 2001 From: Thanhthanh Nguyen Date: Sun, 7 Sep 2025 16:33:17 -0700 Subject: [PATCH 8/9] trying to make it parallalizable on gpu --- .../calibration/calibrate_all_sites.jl | 21 +- experiments/calibration/test_calibration.jl | 757 +++++++++--------- 2 files changed, 395 insertions(+), 383 deletions(-) diff --git a/experiments/calibration/calibrate_all_sites.jl b/experiments/calibration/calibrate_all_sites.jl index 5612dd5c76..897efa0257 100644 --- a/experiments/calibration/calibrate_all_sites.jl +++ b/experiments/calibration/calibrate_all_sites.jl @@ -2,8 +2,14 @@ using ClimaLand using ClimaLand.Artifacts using DelimitedFiles using Logging +using Base.Threads + +using CUDA +using ClimaComms +ClimaComms.@import_required_backends() script = joinpath(pkgdir(ClimaLand), "experiments/calibration/test_calibration.jl") +include(script) # to load dependencies # Get all site IDs from the metadata file fluxnet2015_data_path = ClimaLand.Artifacts.fluxnet2015_data_path() @@ -19,16 +25,17 @@ site_id_col = findfirst(==("site_id"), header) site_ids = table[:, site_id_col] # vector of all site IDs unique_site_ids = unique(site_ids) - - # Run the script for each unique site ID -for site_ID in unique_site_ids - println("Running $script with file: $site_ID") +Threads.@threads for i in 1:length(unique_site_ids) + site_ID = unique_site_ids[i] + + @info "Running $script with file: $site_ID." try - run(`julia --project=.buildkite $script $site_ID`) - println("Success on $site_ID") + # run(`julia --project=.buildkite $script $site_ID`) + calibrate_at_site(site_ID) + @info "Success on $site_ID" catch e - println("Failed on $site_ID: $e") + @info "Failed on $site_ID: $e" continue end end \ No newline at end of file diff --git a/experiments/calibration/test_calibration.jl b/experiments/calibration/test_calibration.jl index 6b776efbeb..192f9f649e 100644 --- a/experiments/calibration/test_calibration.jl +++ b/experiments/calibration/test_calibration.jl @@ -52,407 +52,412 @@ using ClimaAnalysis, GeoMakie, Printf, StatsBase using CSV using DataFrames using DelimitedFiles + + + +rng_seed = 1234 +rng = Random.MersenneTwister(rng_seed) +const FT = Float64 + # ## Configuration and Site Setup # # Configure the experiment parameters and set up the FLUXNET site (US-MOz) with # its specific location, time settings, and atmospheric conditions. # Set random seed for reproducibility and floating point precision -rng_seed = 1234 -rng = Random.MersenneTwister(rng_seed) -const FT = Float64 +function calibrate_at_site(site_ID) + # Initialize land parameters and site configuration. + earth_param_set = LP.LandParameters(FT) + # site_ID = ARGS[1] # should already be defined when including this script + # site_ID_val = FluxnetSimulations.replace_hyphen(site_ID); + + # Get site-specific information: location coordinates, time offset, and sensor + # height. + (; time_offset, lat, long, atmos_h) = + FluxnetSimulations.get_location(site_ID) + + # Get maximum simulation start and end dates in UTC; these must be included in the forcing data range + start_date = DateTime(2011, 3, 1) + stop_date = DateTime(2011, 7, 1) + Δt = 450.0; # seconds + + # ## Domain and Forcing Setup + # + # Create the computational domain and load the necessary forcing data for the + # land surface model. + # ClimaLand includes the domain, forcing, and LAI as + # part of the model, but here we will need to recreate a model + # many times (for each parameter value tried), while the domain, forcing, + # and LAI are held fixed. Because of that, we define them here, once, outside of + # the function that is called to make the model. + + # Create a column domain representing a 2-meter deep soil column with 10 + # vertical layers. + zmin = FT(-2) # 2m depth + zmax = FT(0) # surface + domain = Column(; zlim = (zmin, zmax), nelements = 20, longlat = (long, lat)); + + # Load prescribed atmospheric and radiative forcing from FLUXNET data + forcing = FluxnetSimulations.prescribed_forcing_fluxnet( + site_ID, + lat, + long, + time_offset, + atmos_h, + start_date, + earth_param_set, + FT, + ); + + # Get Leaf Area Index (LAI) data from MODIS satellite observations. + LAI = + ClimaLand.prescribed_lai_modis(domain.space.surface, start_date, stop_date); + + # ## Model Setup + # + # Create an integrated land model that couples canopy, snow, soil, and soil CO2 + # components. This comprehensive model allows us to simulate the full land + # surface system and its interactions. + function model(Vcmax25, g1) + Vcmax25 = FT(Vcmax25) + g1 = FT(g1) + + #md # Set up models; note: we are not using the default soil, which relies on global + #md # maps of parameters to estimate the parameters at the site. + #md # Instead we use parameter more tailored to this site. + prognostic_land_components = (:canopy, :snow, :soil, :soilco2) + #md # Create soil model + retention_parameters = (; + ν = FT(0.55), + θ_r = FT(0.04), + K_sat = FT(4e-7), + hydrology_cm = ClimaLand.Soil.vanGenuchten{FT}(; + α = FT(0.05), + n = FT(2.0), + ), + ) + composition_parameters = + (; ν_ss_om = FT(0.1), ν_ss_quartz = FT(0.1), ν_ss_gravel = FT(0.0)) + soil = ClimaLand.Soil.EnergyHydrology{FT}( + domain, + forcing, + earth_param_set; + prognostic_land_components, + additional_sources = (ClimaLand.RootExtraction{FT}(),), + runoff = ClimaLand.Soil.Runoff.SurfaceRunoff(), + retention_parameters, + composition_parameters, + S_s = FT(1e-2), + ) + surface_domain = ClimaLand.Domains.obtain_surface_domain(domain) + canopy_forcing = (; + atmos = forcing.atmos, + radiation = forcing.radiation, + ground = ClimaLand.PrognosticGroundConditions{FT}(), + ) + #md # Set up photosynthesis using the Farquhar model + photosyn_defaults = + Canopy.clm_photosynthesis_parameters(surface_domain.space.surface) + photosynthesis = Canopy.FarquharModel{FT}( + surface_domain; + photosynthesis_parameters = (; + is_c3 = photosyn_defaults.is_c3, + Vcmax25, + ), + ) + #md # Set up stomatal conductance using the Medlyn model + conductance = Canopy.MedlynConductanceModel{FT}(surface_domain; g1) + + #md # Create canopy model + canopy = Canopy.CanopyModel{FT}( + surface_domain, + canopy_forcing, + LAI, + earth_param_set; + photosynthesis, + conductance, + prognostic_land_components, + ) -# Initialize land parameters and site configuration. -earth_param_set = LP.LandParameters(FT) -site_ID = ARGS[1] -# site_ID_val = FluxnetSimulations.replace_hyphen(site_ID); + #md # Create integrated land model + land_model = + LandModel{FT}(forcing, LAI, earth_param_set, domain, Δt; soil, canopy) -# Get site-specific information: location coordinates, time offset, and sensor -# height. -(; time_offset, lat, long, atmos_h) = - FluxnetSimulations.get_location(site_ID) + #md # Set initial conditions from FLUXNET data + set_ic! = FluxnetSimulations.make_set_fluxnet_initial_conditions( + site_ID, + start_date, + time_offset, + land_model, + ) -# Get maximum simulation start and end dates in UTC; these must be included in the forcing data range -start_date = DateTime(2011, 3, 1) -stop_date = DateTime(2011, 7, 1) -Δt = 450.0; # seconds + #md # Configure diagnostics to output sensible and latent heat fluxes hourly + output_vars = ["lhf"] + diagnostics = ClimaLand.default_diagnostics( + land_model, + start_date; + output_writer = ClimaDiagnostics.Writers.DictWriter(), + output_vars, + average_period = :hourly, + ) -# ## Domain and Forcing Setup -# -# Create the computational domain and load the necessary forcing data for the -# land surface model. -# ClimaLand includes the domain, forcing, and LAI as -# part of the model, but here we will need to recreate a model -# many times (for each parameter value tried), while the domain, forcing, -# and LAI are held fixed. Because of that, we define them here, once, outside of -# the function that is called to make the model. - -# Create a column domain representing a 2-meter deep soil column with 10 -# vertical layers. -zmin = FT(-2) # 2m depth -zmax = FT(0) # surface -domain = Column(; zlim = (zmin, zmax), nelements = 20, longlat = (long, lat)); - -# Load prescribed atmospheric and radiative forcing from FLUXNET data -forcing = FluxnetSimulations.prescribed_forcing_fluxnet( - site_ID, - lat, - long, - time_offset, - atmos_h, - start_date, - earth_param_set, - FT, -); - -# Get Leaf Area Index (LAI) data from MODIS satellite observations. -LAI = - ClimaLand.prescribed_lai_modis(domain.space.surface, start_date, stop_date); - -# ## Model Setup -# -# Create an integrated land model that couples canopy, snow, soil, and soil CO2 -# components. This comprehensive model allows us to simulate the full land -# surface system and its interactions. -function model(Vcmax25, g1) - Vcmax25 = FT(Vcmax25) - g1 = FT(g1) - - #md # Set up models; note: we are not using the default soil, which relies on global - #md # maps of parameters to estimate the parameters at the site. - #md # Instead we use parameter more tailored to this site. - prognostic_land_components = (:canopy, :snow, :soil, :soilco2) - #md # Create soil model - retention_parameters = (; - ν = FT(0.55), - θ_r = FT(0.04), - K_sat = FT(4e-7), - hydrology_cm = ClimaLand.Soil.vanGenuchten{FT}(; - α = FT(0.05), - n = FT(2.0), - ), - ) - composition_parameters = - (; ν_ss_om = FT(0.1), ν_ss_quartz = FT(0.1), ν_ss_gravel = FT(0.0)) - soil = ClimaLand.Soil.EnergyHydrology{FT}( - domain, - forcing, - earth_param_set; - prognostic_land_components, - additional_sources = (ClimaLand.RootExtraction{FT}(),), - runoff = ClimaLand.Soil.Runoff.SurfaceRunoff(), - retention_parameters, - composition_parameters, - S_s = FT(1e-2), - ) - surface_domain = ClimaLand.Domains.obtain_surface_domain(domain) - canopy_forcing = (; - atmos = forcing.atmos, - radiation = forcing.radiation, - ground = ClimaLand.PrognosticGroundConditions{FT}(), - ) - #md # Set up photosynthesis using the Farquhar model - photosyn_defaults = - Canopy.clm_photosynthesis_parameters(surface_domain.space.surface) - photosynthesis = Canopy.FarquharModel{FT}( - surface_domain; - photosynthesis_parameters = (; - is_c3 = photosyn_defaults.is_c3, - Vcmax25, + #md # Create and run the simulation + data_dt = Second(FluxnetSimulations.get_data_dt(site_ID)) + updateat = Array(start_date:Second(Δt):stop_date) + #md # Create and run the simulation + simulation = Simulations.LandSimulation( + start_date, + stop_date, + Δt, + land_model; + set_ic!, + updateat, + user_callbacks = (), + diagnostics, + ) + solve!(simulation) + return simulation + end; + + # ## Observation and Helper Functions + # + # Define the observation function `G` that maps from parameter space to + # observation space, along with supporting functions for data processing: + + # This function runs the model and computes diurnal average of latent heat flux + function G(Vcmax25, g1) + simulation = model(Vcmax25, g1) + lhf = get_lhf(simulation) + observation = + Float64.( + get_diurnal_average( + lhf, + simulation.start_date, + simulation.start_date + Day(20), + stop_date, + ) + ) + return observation + end; + + # Helper function: Extract latent heat flux from simulation diagnostics + function get_lhf(simulation) + return ClimaLand.Diagnostics.diagnostic_as_vectors( + simulation.diagnostics[1].output_writer, + "lhf_1h_average", + ) + end; + + # Helper function: Compute diurnal average of a variable + function get_diurnal_average(var, start_date, spinup_date, stop_date) + (times, data) = var + model_dates = if times isa Vector{DateTime} + times + else + date.(times) + end + spinup_idx = findfirst(spinup_date .<= model_dates) + stop_idx = findlast(model_dates .< stop_date) + model_dates = model_dates[spinup_idx:stop_idx] + data = data[spinup_idx:stop_idx] + + hour_of_day = Hour.(model_dates) + mean_by_hour = [mean(data[hour_of_day .== Hour(i)]) for i in 0:23] + return mean_by_hour + end; + + # ## Experiment Setup + # + # We obtain observations from the FLUXNET site. The dataset contains multiple + # variables, but we will just use latent heat flux in this calibration. + dataset = FluxnetSimulations.get_comparison_data(site_ID, time_offset) + observations = get_diurnal_average( + (dataset.UTC_datetime, dataset.lhf), + start_date, + start_date + Day(20), + stop_date, + ); + + # Define observation noise covariance for the ensemble Kalman process. A flat + # covariance matrix is used here for simplicity. + noise_covariance = 0.05 * EKP.I; + + # ## Prior Distribution and Calibration Configuration + # + # Set up the prior distribution for the parameter and configure the ensemble + # Kalman inversion: + + # Constrained Gaussian prior for Vcmax25 with bounds [0, 1e-3], and g1 with + # bounds [0, 1000]. The first two arguments of each prior are the mean and standard + # deviation of the Gaussian function. + + priors = [ + PD.constrained_gaussian("Vcmax25", 1e-4, 5e-5, 0, 1e-3), + PD.constrained_gaussian("g1", 150, 90, 0, 1000), + ] + prior = PD.combine_distributions(priors); + + # Set the ensemble size and number of iterations + ensemble_size = 10 + N_iterations = 4; + + # ## Ensemble Kalman Inversion + # + # Initialize and run the ensemble Kalman process: + + # Sample the initial parameter ensemble from the prior distribution + initial_ensemble = EKP.construct_initial_ensemble(rng, prior, ensemble_size) + + ensemble_kalman_process = EKP.EnsembleKalmanProcess( + initial_ensemble, + observations, + noise_covariance, + EKP.Inversion(); + scheduler = EKP.DataMisfitController( + terminate_at = Inf, + on_terminate = "continue", ), - ) - #md # Set up stomatal conductance using the Medlyn model - conductance = Canopy.MedlynConductanceModel{FT}(surface_domain; g1) - - #md # Create canopy model - canopy = Canopy.CanopyModel{FT}( - surface_domain, - canopy_forcing, - LAI, - earth_param_set; - photosynthesis, - conductance, - prognostic_land_components, + rng, + ); + + # Run the ensemble of forward models to iteratively update the parameter + # ensemble + Logging.with_logger(SimpleLogger(devnull, Logging.Error)) do + for i in 1:N_iterations + println("Iteration $i") + params_i = EKP.get_ϕ_final(prior, ensemble_kalman_process) + G_ens = hcat([G(params_i[:, j]...) for j in 1:ensemble_size]...) + EKP.update_ensemble!(ensemble_kalman_process, G_ens) + end + end; + + # ## Results Analysis and Visualization + # Get the mean of the final parameter ensemble: + final_mean_params = EKP.get_ϕ_mean_final(prior, ensemble_kalman_process); + + ϕ = EKP.get_ϕ_final(prior, ensemble_kalman_process); + y = EKP.get_obs(ensemble_kalman_process); + metrics = EKP.get_error_metrics(ensemble_kalman_process); + + # Now, let's analyze the calibration results by examining parameter evolution + # and comparing model outputs across iterations. + # + savedir = + joinpath(pkgdir(ClimaLand), "experiments/calibration/out/$(site_ID)/out") + mkpath(savedir) + + # Save the final parameter ensemble and error metrics to CSV files for further analysis + # or record-keeping. + writedlm( + joinpath(savedir, "final_mean_parameters_$(site_ID).csv"), + final_mean_params, + ',', ) - #md # Create integrated land model - land_model = - LandModel{FT}(forcing, LAI, earth_param_set, domain, Δt; soil, canopy) - - #md # Set initial conditions from FLUXNET data - set_ic! = FluxnetSimulations.make_set_fluxnet_initial_conditions( - site_ID, - start_date, - time_offset, - land_model, + writedlm( + joinpath(savedir, "final_parameter_ensemble_$(site_ID).csv"), + ϕ, + ',', ) - #md # Configure diagnostics to output sensible and latent heat fluxes hourly - output_vars = ["lhf"] - diagnostics = ClimaLand.default_diagnostics( - land_model, - start_date; - output_writer = ClimaDiagnostics.Writers.DictWriter(), - output_vars, - average_period = :hourly, + writedlm( + joinpath(savedir, "latest_data_vector_$(site_ID).csv"), + y, + ',', ) - #md # Create and run the simulation - data_dt = Second(FluxnetSimulations.get_data_dt(site_ID)) - updateat = Array(start_date:Second(Δt):stop_date) - #md # Create and run the simulation - simulation = Simulations.LandSimulation( - start_date, - stop_date, - Δt, - land_model; - set_ic!, - updateat, - user_callbacks = (), - diagnostics, - ) - solve!(simulation) - return simulation -end; + df = DataFrame(metrics) + CSV.write(savedir, "error_metrics_$(site_ID).csv", df) -# ## Observation and Helper Functions -# -# Define the observation function `G` that maps from parameter space to -# observation space, along with supporting functions for data processing: - -# This function runs the model and computes diurnal average of latent heat flux -function G(Vcmax25, g1) - simulation = model(Vcmax25, g1) - lhf = get_lhf(simulation) - observation = - Float64.( - get_diurnal_average( - lhf, - simulation.start_date, - simulation.start_date + Day(20), - stop_date, - ) - ) - return observation -end; - -# Helper function: Extract latent heat flux from simulation diagnostics -function get_lhf(simulation) - return ClimaLand.Diagnostics.diagnostic_as_vectors( - simulation.diagnostics[1].output_writer, - "lhf_1h_average", - ) -end; - -# Helper function: Compute diurnal average of a variable -function get_diurnal_average(var, start_date, spinup_date, stop_date) - (times, data) = var - model_dates = if times isa Vector{DateTime} - times - else - date.(times) - end - spinup_idx = findfirst(spinup_date .<= model_dates) - stop_idx = findlast(model_dates .< stop_date) - model_dates = model_dates[spinup_idx:stop_idx] - data = data[spinup_idx:stop_idx] + # Plot the parameter ensemble evolution over iterations to visualize + # convergence: - hour_of_day = Hour.(model_dates) - mean_by_hour = [mean(data[hour_of_day .== Hour(i)]) for i in 0:23] - return mean_by_hour -end; + dim_size = sum(length.(EKP.batch(prior))) + fig = CairoMakie.Figure(size = ((dim_size + 1) * 500, 500)) -# ## Experiment Setup -# -# We obtain observations from the FLUXNET site. The dataset contains multiple -# variables, but we will just use latent heat flux in this calibration. -dataset = FluxnetSimulations.get_comparison_data(site_ID, time_offset) -observations = get_diurnal_average( - (dataset.UTC_datetime, dataset.lhf), - start_date, - start_date + Day(20), - stop_date, -); - -# Define observation noise covariance for the ensemble Kalman process. A flat -# covariance matrix is used here for simplicity. -noise_covariance = 0.05 * EKP.I; - -# ## Prior Distribution and Calibration Configuration -# -# Set up the prior distribution for the parameter and configure the ensemble -# Kalman inversion: + for i in 1:dim_size + EKP.Visualize.plot_ϕ_over_iters( + fig[1, i], + ensemble_kalman_process, + prior, + i, + ) + end -# Constrained Gaussian prior for Vcmax25 with bounds [0, 1e-3], and g1 with -# bounds [0, 1000]. The first two arguments of each prior are the mean and standard -# deviation of the Gaussian function. + # TODO save numerical error over iterations + EKP.Visualize.plot_error_over_iters( + fig[1, dim_size + 1], + ensemble_kalman_process, + ) + CairoMakie.save("$(savedir)/constrained_params_and_error.png", fig) + # fig + + # Compare the model output between the first and last iterations to assess + # improvement: + fig = CairoMakie.Figure(size = (900, 400)) + + first_G_ensemble = EKP.get_g(ensemble_kalman_process, 1) + last_iter = EKP.get_N_iterations(ensemble_kalman_process) + last_G_ensemble = EKP.get_g(ensemble_kalman_process, last_iter) + n_ens = EKP.get_N_ens(ensemble_kalman_process) + + writedlm( + joinpath(savedir, "first_iteration_G_ensemble_$(site_ID).csv"), + first_G_ensemble, + ',', + ) -priors = [ - PD.constrained_gaussian("Vcmax25", 1e-4, 5e-5, 0, 1e-3), - PD.constrained_gaussian("g1", 150, 90, 0, 1000), -] -prior = PD.combine_distributions(priors); + writedlm( + joinpath(savedir, "last_iteration_G_ensemble_$(site_ID).csv"), + last_G_ensemble, + ',', + ) -# Set the ensemble size and number of iterations -ensemble_size = 10 -N_iterations = 4; + ax = Axis( + fig[1, 1]; + title = "G ensemble: first vs last iteration (n = $(n_ens), iters 1 vs $(last_iter))", + xlabel = "Observation index", + ylabel = "G", + ) -# ## Ensemble Kalman Inversion -# -# Initialize and run the ensemble Kalman process: - -# Sample the initial parameter ensemble from the prior distribution -initial_ensemble = EKP.construct_initial_ensemble(rng, prior, ensemble_size) - -ensemble_kalman_process = EKP.EnsembleKalmanProcess( - initial_ensemble, - observations, - noise_covariance, - EKP.Inversion(); - scheduler = EKP.DataMisfitController( - terminate_at = Inf, - on_terminate = "continue", - ), - rng, -); - -# Run the ensemble of forward models to iteratively update the parameter -# ensemble -Logging.with_logger(SimpleLogger(devnull, Logging.Error)) do - for i in 1:N_iterations - println("Iteration $i") - params_i = EKP.get_ϕ_final(prior, ensemble_kalman_process) - G_ens = hcat([G(params_i[:, j]...) for j in 1:ensemble_size]...) - EKP.update_ensemble!(ensemble_kalman_process, G_ens) + # Plot model output of first vs last iteration ensemble + for g in eachcol(first_G_ensemble) + lines!(ax, 1:length(g), g; color = (:red, 0.6), linewidth = 1.5) end -end; -# ## Results Analysis and Visualization -# Get the mean of the final parameter ensemble: -final_mean_params = EKP.get_ϕ_mean_final(prior, ensemble_kalman_process); - -ϕ = EKP.get_ϕ_final(prior, ensemble_kalman_process); -y = EKP.get_obs(ensemble_kalman_process); -metrics = EKP.get_error_metrics(ensemble_kalman_process); + for g in eachcol(last_G_ensemble) + lines!(ax, 1:length(g), g; color = (:blue, 0.6), linewidth = 1.5) + end -# Now, let's analyze the calibration results by examining parameter evolution -# and comparing model outputs across iterations. -# -savedir = - joinpath(pkgdir(ClimaLand), "experiments/calibration/out/$(site_ID)/out") -mkpath(savedir) - -# Save the final parameter ensemble and error metrics to CSV files for further analysis -# or record-keeping. -writedlm( - joinpath(savedir, "final_mean_parameters_$(site_ID).csv"), - final_mean_params, - ',', -) - -writedlm( - joinpath(savedir, "final_parameter_ensemble_$(site_ID).csv"), - ϕ, - ',', -) - -writedlm( - joinpath(savedir, "latest_data_vector_$(site_ID).csv"), - y, - ',', -) - -df = DataFrame(metrics) -CSV.write(savedir, "error_metrics_$(site_ID).csv", df) - -# Plot the parameter ensemble evolution over iterations to visualize -# convergence: - -dim_size = sum(length.(EKP.batch(prior))) -fig = CairoMakie.Figure(size = ((dim_size + 1) * 500, 500)) - -for i in 1:dim_size - EKP.Visualize.plot_ϕ_over_iters( - fig[1, i], - ensemble_kalman_process, - prior, - i, + # --- Add shaded noise around observations --- + xvals = 1:length(observations) + + # Say you have obs variance as `obs_var` (same length as `observations`) + # If you don’t, but only have noise covariance = 0.05*I, you can just use sqrt(0.05) + obs_std = sqrt.(0.05) .* length(observations) # or your own per-point std + + band!(ax, xvals, + observations .- 2 .* obs_std, + observations .+ 2 .* obs_std; + color = (:gray, 0.3), + transparency = true) + + # now lines on top + lines!(ax, xvals, observations; color = (:black, 0.8), linewidth = 3) + + axislegend( + ax, + [ + LineElement(color = :red, linewidth = 2), + LineElement(color = :blue, linewidth = 2), + LineElement(color = :black, linewidth = 4), + ], + ["First ensemble", "Last ensemble", "Observations"]; + position = :rb, + framevisible = false, ) -end - -# TODO save numerical error over iterations -EKP.Visualize.plot_error_over_iters( - fig[1, dim_size + 1], - ensemble_kalman_process, -) -CairoMakie.save("$(savedir)/constrained_params_and_error.png", fig) -# fig - -# Compare the model output between the first and last iterations to assess -# improvement: -fig = CairoMakie.Figure(size = (900, 400)) - -first_G_ensemble = EKP.get_g(ensemble_kalman_process, 1) -last_iter = EKP.get_N_iterations(ensemble_kalman_process) -last_G_ensemble = EKP.get_g(ensemble_kalman_process, last_iter) -n_ens = EKP.get_N_ens(ensemble_kalman_process) - -writedlm( - joinpath(savedir, "first_iteration_G_ensemble_$(site_ID).csv"), - first_G_ensemble, - ',', -) - -writedlm( - joinpath(savedir, "last_iteration_G_ensemble_$(site_ID).csv"), - last_G_ensemble, - ',', -) - -ax = Axis( - fig[1, 1]; - title = "G ensemble: first vs last iteration (n = $(n_ens), iters 1 vs $(last_iter))", - xlabel = "Observation index", - ylabel = "G", -) - -# Plot model output of first vs last iteration ensemble -for g in eachcol(first_G_ensemble) - lines!(ax, 1:length(g), g; color = (:red, 0.6), linewidth = 1.5) -end - -for g in eachcol(last_G_ensemble) - lines!(ax, 1:length(g), g; color = (:blue, 0.6), linewidth = 1.5) -end - -# --- Add shaded noise around observations --- -xvals = 1:length(observations) - -# Say you have obs variance as `obs_var` (same length as `observations`) -# If you don’t, but only have noise covariance = 0.05*I, you can just use sqrt(0.05) -obs_std = sqrt.(0.05) .* length(observations) # or your own per-point std - -band!(ax, xvals, - observations .- 2 .* obs_std, - observations .+ 2 .* obs_std; - color = (:gray, 0.3), - transparency = true) - -# now lines on top -lines!(ax, xvals, observations; color = (:black, 0.8), linewidth = 3) - -axislegend( - ax, - [ - LineElement(color = :red, linewidth = 2), - LineElement(color = :blue, linewidth = 2), - LineElement(color = :black, linewidth = 4), - ], - ["First ensemble", "Last ensemble", "Observations"]; - position = :rb, - framevisible = false, -) - -CairoMakie.resize_to_layout!(fig) -CairoMakie.save("$(savedir)/G_first_and_last.png", fig) -# fig \ No newline at end of file + + CairoMakie.resize_to_layout!(fig) + CairoMakie.save("$(savedir)/G_first_and_last.png", fig) + # fig +end \ No newline at end of file From 7351f814520d509b75e2eafa41cbcf8ec6dda71e Mon Sep 17 00:00:00 2001 From: ThanhNguyen428 Date: Tue, 28 Oct 2025 22:03:19 -0700 Subject: [PATCH 9/9] cleaning up for a branch clone --- .buildkite/Manifest-v1.11.toml | 4 +- .buildkite/Project.toml | 3 +- Artifacts.toml | 191 -------- .../calibration/calibrate_all_sites.jl | 31 -- .../calibration/calibrate_all_sites.jl | 41 -- experiments/calibration/test_calibration.jl | 463 ------------------ experiments/integrated/fluxnet/run_fluxnet.jl | 17 +- ext/fluxnet_simulations/generic_site.jl | 22 +- .../get_fluxnet_metadata.jl | 6 +- ext/land_sim_vis/plotting_utils.jl | 2 +- src/Artifacts.jl | 17 +- test/integrated/fluxnet_sim.jl | 6 +- 12 files changed, 44 insertions(+), 759 deletions(-) delete mode 100644 Artifacts.toml delete mode 100644 docs/src/tutorials/calibration/calibrate_all_sites.jl delete mode 100644 experiments/calibration/calibrate_all_sites.jl delete mode 100644 experiments/calibration/test_calibration.jl diff --git a/.buildkite/Manifest-v1.11.toml b/.buildkite/Manifest-v1.11.toml index 0b0171984f..c990df681f 100644 --- a/.buildkite/Manifest-v1.11.toml +++ b/.buildkite/Manifest-v1.11.toml @@ -1,8 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.11.6" +julia_version = "1.11.7" manifest_format = "2.0" -project_hash = "0257c8e4fc6fd73b316df90fb308c7d86009613a" +project_hash = "53802e4fba14dac3d5e5fb1149db2ae8e60ebb4d" [[deps.ADTypes]] git-tree-sha1 = "7927b9af540ee964cc5d1b73293f1eb0b761a3a1" diff --git a/.buildkite/Project.toml b/.buildkite/Project.toml index 80f610aad9..8d9a718434 100644 --- a/.buildkite/Project.toml +++ b/.buildkite/Project.toml @@ -37,6 +37,7 @@ ProfileCanvas = "efd6af41-a80b-495e-886c-e51b0c7d77a3" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" @@ -52,5 +53,5 @@ ClimaDiagnostics = "0.2.13" ClimaTimeSteppers = "0.7, 0.8" EnsembleKalmanProcesses = "2.5.0" Flux = "0.15" -GeoMakie = "< 0.7.13" # v0.7.13 causes infinite recursion with GridPositions +GeoMakie = "< 0.7.13" Statistics = "1" diff --git a/Artifacts.toml b/Artifacts.toml deleted file mode 100644 index 8e79b8b354..0000000000 --- a/Artifacts.toml +++ /dev/null @@ -1,191 +0,0 @@ -[landsea_mask_30arcseconds] -git-tree-sha1 = "517c925535981f9d17ac9e7a65e2c35c3ce9597d" - - [[landsea_mask_30arcseconds.download]] - sha256 = "11eaf6b7606056d198f3735f7b76ead76269cced186de78075f415101cf39088" - url = "https://caltech.box.com/shared/static/z013ay1k6oqofnw69bv5held9x374sux.gz" - -[landsea_mask_60arcseconds] -git-tree-sha1 = "10ed7ec61def091c44545825cfb4d9599216c3c8" - - [[landsea_mask_60arcseconds.download]] - sha256 = "eb88455a39a00bfadb77effbc31a4df561158c0ca1578476edb683d45e3c49da" - url = "https://caltech.box.com/shared/static/vivf36pih7e0ttahqarc2nd5jxb4q8jr.gz" - -[landsea_mask_1deg] -git-tree-sha1 = "cbbc6b3752d9cb9b667ec33cfbeb46819f8db418" - - [[landsea_mask_1deg.download]] - sha256 = "3722b553c2fdf28a6574aea2e0b167d16ab050f34e5969ada45625b3a3ecb6da" - url = "https://caltech.box.com/shared/static/b3u4dv0dsoswvqgp8y63bzj7awbhwztd.gz" - -[soil_ic_2008_50m] -git-tree-sha1 = "7b0f0572f4620a8ba2a7ae816a547d2dcfb9eb9b" - - [[soil_ic_2008_50m.download]] - sha256 = "b0080d605585072b88871bc8987f615c9dfc35e64ad16e00a52087d879fc5c3a" - url = "https://caltech.box.com/shared/static/atiztw8uu2i78g2vd02n8onb68cer9tn.gz" - -[soilgrids] -git-tree-sha1 = "f919da748cc973519fb136f4225cf343f5c3e41f" -[soilgrids_lowres] -git-tree-sha1 = "6ce1717bb8b7aa767380f1053b9d8d3a4e4edaa6" - - [[soilgrids_lowres.download]] - sha256 = "27d49fce7f4a5e2b81bbd29e8a7caf08f568a2df53e3af817126363decfd5a7c" - url = "https://caltech.box.com/shared/static/99aw2gce2k65bdu0h8jkfxq4vin08gi2.gz" - -[sw_albedo] -git-tree-sha1 = "136f1db3ed969614fb589c5250545a3ed1e8aaab" - - [[sw_albedo.download]] - sha256 = "777371543cc3285c0779e7ecd925fe782b5b45e0ed27f82260827ac8ebe6a700" - url = "https://caltech.box.com/shared/static/of1admpndmikoumtgk5j3yvt92v71awk.gz" - -["clm_data_0.9x1.25"] -git-tree-sha1 = "38488be5dd476890deb6fbd2b26e3e1cee2d449e" - - [["clm_data_0.9x1.25".download]] - sha256 = "cb5ead78868243879261d26cbb42d665bfd98135ebbebc052487d4d13ebf3a98" - url = "https://caltech.box.com/shared/static/6ms2qzyiekorksewmfyie6mbonlq5ctp.gz" -["clm_data_0.125x0.125"] -git-tree-sha1 = "b8816f0af2c6d182111a156c3b8eed2ed83efea7" - - [["clm_data_0.125x0.125".download]] - sha256 = "993b5b97b6ec22ec484db4e523b65a231d736ead8b5e0759ad3abd00050bb58c" - url = "https://caltech.box.com/shared/static/fx9138ysguyg4g65l4tvibbk16er0nqo.gz" - -[era5_land_forcing_data2008] -git-tree-sha1 = "93d2e93f491e77cb8fba2a1b8b3946f38bde469e" -[era5_land_forcing_data2008_lowres] -git-tree-sha1 = "fd826542e6ab1758c8d5a2cc6dbf1f67b47aadee" - - [[era5_land_forcing_data2008_lowres.download]] - sha256 = "c5ea5acfcab2e69c26d6ef8ba20e9e6dc018e5baa9c48d0df69cf0a5f9bf3399" - url = "https://caltech.box.com/shared/static/6zk5fa3ka5ib62jsrs6tx35sqwe8xv2k.gz" -[era5_land_forcing_data2008_lai] -git-tree-sha1 = "df5ec5af1ea0793c0e7d49e60db5a938098e174f" - - [[era5_land_forcing_data2008_lai.download]] - sha256 = "a83bf28a7a864db285e41289031a9761c2a426373a6efecd04b866a0b5374cef" - url = "https://caltech.box.com/shared/static/m82r23e67x04a3hb6p79dfltlumvg9z0.gz" - -[modis_lai] -git-tree-sha1 = "81d4bc1b22e94d66eec3d80a2d21b5dd5303bf8f" - - [[modis_lai.download]] - sha256 = "4a244428fe480aa05ff863af69299d2ed8ca66ed324b9e1bffcb440792fb4cce" - url = "https://caltech.box.com/shared/static/99ktg026kciezykyegsawd1gjwj1du74.gz" - -[land_albedo] -git-tree-sha1 = "3f07b70ab43a05123e93753a45e094dcf7a66b5b" - -[snowmip] -git-tree-sha1 = "230f00823299e7e33d276d9b57cba985fc04eaae" - - [[snowmip.download]] - sha256 = "598b525c7421fabf06de4c8d4d7e05d50b7655ad2387d98c4335e90fd5633d61" - url = "https://caltech.box.com/shared/static/0vqhet1hsj36sgmy0ajrcrb24cihcfd7.gz" - -[soil_params_Gupta2020_2022] -git-tree-sha1 = "c10e8f877603bfa2c018eb2dd0c51dd078abe5df" - - [[soil_params_Gupta2020_2022.download]] - sha256 = "6205d96ae7724c469eef5cd1cf3b99e186637427b6fb65ecc659fb808fc6cacf" - url = "https://caltech.box.com/shared/static/7c6yx62tzjivxfmhqyy03r6e9l1glyu1.gz" - -[modis_clumping_index] -git-tree-sha1 = "b849eb95c09190095e7bf021494ddeda8858af01" - - [[modis_clumping_index.download]] - sha256 = "e4c766a93a673e5dc22540687ef5616416d65bb13a0f4a67789b95d49ccbb158" - url = "https://caltech.box.com/shared/static/ec2y3k5kqpl9wjvtsx3584wgkp5q8dyw.gz" - -[bonan_richards_eqn] -git-tree-sha1 = "043f9354b961fd3ef6ac4cf71d0c99930e177a92" - [[bonan_richards_eqn.download]] - sha256 = "50f1739dfd8193742488f249bd1ba8a157dfc7c1d8a27392ba7edbe9b4c0db03" - url = "https://caltech.box.com/shared/static/4jcz9c1lp6rt750q28kx3qvihi8393tv.gz" - -[fluxnet_sites] -git-tree-sha1 = "2fc70601badf6f83dee2b84ba9c386ad041de8e2" - - [[fluxnet_sites.download]] - sha256 = "f05b4c01b57afe9c0f59095b39cea1c0cf46b20deecd5c7afb44f474d9cd9966" - url = "https://caltech.box.com/shared/static/otrr2y0rgjct7hqhmq214nb8qjsvqj5p.gz" - -[topmodel] -git-tree-sha1 = "10855ff977cc948a92c7a83915895d27650649a1" - - [[topmodel.download]] - sha256 = "03de8f42fedf3836c6792571f04b55fbaca321919ecd761b846286dc56bfa2e8" - url = "https://caltech.box.com/shared/static/bfjxp72vre6b5ncm7jqf9w8hyj4ab4yr.gz" - -[lehmann2008_evaporation] -git-tree-sha1 = "49386718eda19d2fa00ae6981d83948c6664057a" - - [[lehmann2008_evaporation.download]] - sha256 = "e9603add174709a8eb7f70487bdc0f519843c03d83cb40cf07acd428f0d9d881" - url = "https://caltech.box.com/shared/static/khhop0nvuumrtzdakemaab8hsq8emh8e.gz" - -[ilamb_data] -git-tree-sha1 = "839224a62b59d73073bdb9a5c55d3dc75e30fe33" - - [[ilamb_data.download]] - sha256 = "64a9a344ebfbb0113014178a1f93a655c401431565907c07fd33aff8860b62d6" - url = "https://caltech.box.com/shared/static/eii2bfwfp47axfeuysgxlgzbczz27u5g.gz" - -[era5_lai_covers] -git-tree-sha1 = "1e77fc84a027da43f1bf6ebe0a53849595a2be22" - - [[era5_lai_covers.download]] - sha256 = "0331e3a2fc2cb53b0f1d9c5a84f49989da0b506258edc229da161e9ca4a61457" - url = "https://caltech.box.com/shared/static/uzuzk5wh8r9q477jbrusfg1sr0ytzem9.gz" - -[earth_orography_30arcseconds] -git-tree-sha1 = "03dd08fcbf363ed055a176dd7adefb60ff1c3493" - -[earth_orography_60arcseconds] -git-tree-sha1 = "fe19d8dbe7a18ff39588e1f718014b0479d9c0f7" -lazy = true - - [[earth_orography_60arcseconds.download]] - sha256 = "eca66c0701d1c2b9e271742314915ffbf4a0fae92709df611c323f38e019966e" - url = "https://caltech.box.com/shared/static/4asrxcgl6xsgenfcug9p0wkkyhtqilgk.gz" - -[bedrock_depth_30arcseconds] -git-tree-sha1 = "cf5d2f9c2e04f930cd30c6ea1a063078113e2a22" - -[bedrock_depth_60arcseconds] -git-tree-sha1 = "4cc329093b0dbc3443db36e7c01133267ac94d37" - - [[bedrock_depth_60arcseconds.download]] - sha256 = "6eb359e900141de837c41c0082257f7cf082e06e1371d45b07b1078f07e96c53" - url = "https://caltech.box.com/shared/static/k1zq0w1psqu8dk0s9jotrfuu6hm8kxfp.gz" - -[forty_yrs_era5_land_forcing_data] -git-tree-sha1 = "f269a0b057b9f438b4caafdef17da73746310787" - -[era5_monthly_averages_surface_single_level_1979_2024] -git-tree-sha1 = "6cddb07eeee2dd46dc5a19e9b2f706302ddba2c9" - - [[era5_monthly_averages_surface_single_level_1979_2024.download]] - sha256 = "46b422722d98c89c6bc0b8641bff259db1caee253f45389ddf9eb9c2d31ed605" - url = "https://caltech.box.com/shared/static/jbgtyt6oq9lxvk8il5zzck6k581q7f3k.gz" - -[huang_van_genuchten_data] -git-tree-sha1 = "39a344685d74346f89c4858c183565554494fc01" - - [[huang_van_genuchten_data.download]] - sha256 = "1a31d971065bec12df089d95a6f385b44bf12fd136428fa7cb1456da87bd4fef" - url = "https://caltech.box.com/shared/static/fjq0f57qurctx9w5fc6vj8ml66wi8r73.gz" - -[mizoguchi_soil_freezing_data] -git-tree-sha1 = "c35ba0e899040cb8153226ab751f69100f475d39" - - [[mizoguchi_soil_freezing_data.download]] - sha256 = "0027cc080ba45ba33dc790b176ec2854353ce7dce4eae4bef72963b0dd944e0b" - url = "https://caltech.box.com/shared/static/tn1bnqjmegyetw5kzd2ixq5pbnb05s3u.gz" - -[fluxnet2015] -git-tree-sha1 = "94d6eb8e3bc5fc361a43597ab4fb6941bdbe2850" \ No newline at end of file diff --git a/docs/src/tutorials/calibration/calibrate_all_sites.jl b/docs/src/tutorials/calibration/calibrate_all_sites.jl deleted file mode 100644 index 30492fa6c7..0000000000 --- a/docs/src/tutorials/calibration/calibrate_all_sites.jl +++ /dev/null @@ -1,31 +0,0 @@ -using ClimaLand -using ClimaLand.Artifacts -using DelimitedFiles - -# folder = "path/to/folder" -script = "test_calibration.jl" - -fluxnet2015_data_path = ClimaLand.Artifacts.fluxnet2015_data_path() -fluxnet2015_metadata_path = - joinpath(fluxnet2015_data_path, "metadata_DD_clean.csv") - - -data = readdlm(fluxnet2015_metadata_path, ',', header=true) - -header = data[2] # column names -table = data[1] # the actual rows - -site_id_col = findfirst(==("SITE_ID"), header) -site_ids = table[:, site_id_col] # vector of all site IDs -unique_site_ids = unique(site_ids) - -for f in readdir(folder; join=true) - try - println("Running $script with file: $f") - run(`julia $script $f`) - println("✅ Success on $f") - catch e - println("❌ Failed on $f: $e") - continue - end -end \ No newline at end of file diff --git a/experiments/calibration/calibrate_all_sites.jl b/experiments/calibration/calibrate_all_sites.jl deleted file mode 100644 index 897efa0257..0000000000 --- a/experiments/calibration/calibrate_all_sites.jl +++ /dev/null @@ -1,41 +0,0 @@ -using ClimaLand -using ClimaLand.Artifacts -using DelimitedFiles -using Logging -using Base.Threads - -using CUDA -using ClimaComms -ClimaComms.@import_required_backends() - -script = joinpath(pkgdir(ClimaLand), "experiments/calibration/test_calibration.jl") -include(script) # to load dependencies - -# Get all site IDs from the metadata file -fluxnet2015_data_path = ClimaLand.Artifacts.fluxnet2015_data_path() -fluxnet2015_metadata_path = - joinpath(fluxnet2015_data_path, "metadata_DD_clean.csv") - -data = readdlm(fluxnet2015_metadata_path, ',', header=true) - -header = data[2] # column names -table = data[1] # the actual rows - -site_id_col = findfirst(==("site_id"), header) -site_ids = table[:, site_id_col] # vector of all site IDs -unique_site_ids = unique(site_ids) - -# Run the script for each unique site ID -Threads.@threads for i in 1:length(unique_site_ids) - site_ID = unique_site_ids[i] - - @info "Running $script with file: $site_ID." - try - # run(`julia --project=.buildkite $script $site_ID`) - calibrate_at_site(site_ID) - @info "Success on $site_ID" - catch e - @info "Failed on $site_ID: $e" - continue - end -end \ No newline at end of file diff --git a/experiments/calibration/test_calibration.jl b/experiments/calibration/test_calibration.jl deleted file mode 100644 index 192f9f649e..0000000000 --- a/experiments/calibration/test_calibration.jl +++ /dev/null @@ -1,463 +0,0 @@ -# # Model Site-Level Calibration Tutorial Using Observations -# -# In this tutorial we will calibrate the Vcmax25 and g1 parameters using latent heat -# flux observations from the FLUXNET site (US-MOz). -# -# ## Overview -# -# The tutorial covers: -# 1. Setting up a land surface model for a FLUXNET site (US-MOz) -# 2. Obtaining the observation dataset -# 3. Implementing Ensemble Kalman Inversion to calibrate Vcmax25 and g1 -# 4. Analyzing the calibration results -# -#nb # ## Prerequisites -#nb # -#nb # First, ensure you have the required packages installed: - -#nb ENV["JULIA_PKG_PRECOMPILE_AUTO"]=0 -#nb using Pkg -#nb required_pkgs = ["ClimaLand", "ClimaDiagnostics", "CairoMakie", -#nb "EnsembleKalmanProcesses", "Random", "Logging"] -#nb Pkg.add(required_pkgs) -#nb plotting_pkgs = ["ClimaAnalysis", "GeoMakie", "Printf", "StatsBase"] -#nb Pkg.add(plotting_pkgs) -#nb Pkg.instantiate() -#nb Pkg.precompile() - -# ## Setup and Imports -# -# Load all the necessary packages for land surface modeling, diagnostics, -# plotting, and ensemble methods: - -using ClimaLand -using ClimaLand.Domains: Column -using ClimaLand.Canopy -using ClimaLand.Simulations -import ClimaLand.FluxnetSimulations as FluxnetSimulations -import ClimaLand.Parameters as LP -import ClimaLand.LandSimVis as LandSimVis -import ClimaDiagnostics -import ClimaUtilities.TimeManager: date -import EnsembleKalmanProcesses as EKP -import EnsembleKalmanProcesses.ParameterDistributions as PD -using CairoMakie -CairoMakie.activate!() -using Statistics -using Logging -import Random -using Dates -using ClimaAnalysis, GeoMakie, Printf, StatsBase - -using CSV -using DataFrames -using DelimitedFiles - - - -rng_seed = 1234 -rng = Random.MersenneTwister(rng_seed) -const FT = Float64 - -# ## Configuration and Site Setup -# -# Configure the experiment parameters and set up the FLUXNET site (US-MOz) with -# its specific location, time settings, and atmospheric conditions. - -# Set random seed for reproducibility and floating point precision -function calibrate_at_site(site_ID) - # Initialize land parameters and site configuration. - earth_param_set = LP.LandParameters(FT) - # site_ID = ARGS[1] # should already be defined when including this script - # site_ID_val = FluxnetSimulations.replace_hyphen(site_ID); - - # Get site-specific information: location coordinates, time offset, and sensor - # height. - (; time_offset, lat, long, atmos_h) = - FluxnetSimulations.get_location(site_ID) - - # Get maximum simulation start and end dates in UTC; these must be included in the forcing data range - start_date = DateTime(2011, 3, 1) - stop_date = DateTime(2011, 7, 1) - Δt = 450.0; # seconds - - # ## Domain and Forcing Setup - # - # Create the computational domain and load the necessary forcing data for the - # land surface model. - # ClimaLand includes the domain, forcing, and LAI as - # part of the model, but here we will need to recreate a model - # many times (for each parameter value tried), while the domain, forcing, - # and LAI are held fixed. Because of that, we define them here, once, outside of - # the function that is called to make the model. - - # Create a column domain representing a 2-meter deep soil column with 10 - # vertical layers. - zmin = FT(-2) # 2m depth - zmax = FT(0) # surface - domain = Column(; zlim = (zmin, zmax), nelements = 20, longlat = (long, lat)); - - # Load prescribed atmospheric and radiative forcing from FLUXNET data - forcing = FluxnetSimulations.prescribed_forcing_fluxnet( - site_ID, - lat, - long, - time_offset, - atmos_h, - start_date, - earth_param_set, - FT, - ); - - # Get Leaf Area Index (LAI) data from MODIS satellite observations. - LAI = - ClimaLand.prescribed_lai_modis(domain.space.surface, start_date, stop_date); - - # ## Model Setup - # - # Create an integrated land model that couples canopy, snow, soil, and soil CO2 - # components. This comprehensive model allows us to simulate the full land - # surface system and its interactions. - function model(Vcmax25, g1) - Vcmax25 = FT(Vcmax25) - g1 = FT(g1) - - #md # Set up models; note: we are not using the default soil, which relies on global - #md # maps of parameters to estimate the parameters at the site. - #md # Instead we use parameter more tailored to this site. - prognostic_land_components = (:canopy, :snow, :soil, :soilco2) - #md # Create soil model - retention_parameters = (; - ν = FT(0.55), - θ_r = FT(0.04), - K_sat = FT(4e-7), - hydrology_cm = ClimaLand.Soil.vanGenuchten{FT}(; - α = FT(0.05), - n = FT(2.0), - ), - ) - composition_parameters = - (; ν_ss_om = FT(0.1), ν_ss_quartz = FT(0.1), ν_ss_gravel = FT(0.0)) - soil = ClimaLand.Soil.EnergyHydrology{FT}( - domain, - forcing, - earth_param_set; - prognostic_land_components, - additional_sources = (ClimaLand.RootExtraction{FT}(),), - runoff = ClimaLand.Soil.Runoff.SurfaceRunoff(), - retention_parameters, - composition_parameters, - S_s = FT(1e-2), - ) - surface_domain = ClimaLand.Domains.obtain_surface_domain(domain) - canopy_forcing = (; - atmos = forcing.atmos, - radiation = forcing.radiation, - ground = ClimaLand.PrognosticGroundConditions{FT}(), - ) - #md # Set up photosynthesis using the Farquhar model - photosyn_defaults = - Canopy.clm_photosynthesis_parameters(surface_domain.space.surface) - photosynthesis = Canopy.FarquharModel{FT}( - surface_domain; - photosynthesis_parameters = (; - is_c3 = photosyn_defaults.is_c3, - Vcmax25, - ), - ) - #md # Set up stomatal conductance using the Medlyn model - conductance = Canopy.MedlynConductanceModel{FT}(surface_domain; g1) - - #md # Create canopy model - canopy = Canopy.CanopyModel{FT}( - surface_domain, - canopy_forcing, - LAI, - earth_param_set; - photosynthesis, - conductance, - prognostic_land_components, - ) - - #md # Create integrated land model - land_model = - LandModel{FT}(forcing, LAI, earth_param_set, domain, Δt; soil, canopy) - - #md # Set initial conditions from FLUXNET data - set_ic! = FluxnetSimulations.make_set_fluxnet_initial_conditions( - site_ID, - start_date, - time_offset, - land_model, - ) - - #md # Configure diagnostics to output sensible and latent heat fluxes hourly - output_vars = ["lhf"] - diagnostics = ClimaLand.default_diagnostics( - land_model, - start_date; - output_writer = ClimaDiagnostics.Writers.DictWriter(), - output_vars, - average_period = :hourly, - ) - - #md # Create and run the simulation - data_dt = Second(FluxnetSimulations.get_data_dt(site_ID)) - updateat = Array(start_date:Second(Δt):stop_date) - #md # Create and run the simulation - simulation = Simulations.LandSimulation( - start_date, - stop_date, - Δt, - land_model; - set_ic!, - updateat, - user_callbacks = (), - diagnostics, - ) - solve!(simulation) - return simulation - end; - - # ## Observation and Helper Functions - # - # Define the observation function `G` that maps from parameter space to - # observation space, along with supporting functions for data processing: - - # This function runs the model and computes diurnal average of latent heat flux - function G(Vcmax25, g1) - simulation = model(Vcmax25, g1) - lhf = get_lhf(simulation) - observation = - Float64.( - get_diurnal_average( - lhf, - simulation.start_date, - simulation.start_date + Day(20), - stop_date, - ) - ) - return observation - end; - - # Helper function: Extract latent heat flux from simulation diagnostics - function get_lhf(simulation) - return ClimaLand.Diagnostics.diagnostic_as_vectors( - simulation.diagnostics[1].output_writer, - "lhf_1h_average", - ) - end; - - # Helper function: Compute diurnal average of a variable - function get_diurnal_average(var, start_date, spinup_date, stop_date) - (times, data) = var - model_dates = if times isa Vector{DateTime} - times - else - date.(times) - end - spinup_idx = findfirst(spinup_date .<= model_dates) - stop_idx = findlast(model_dates .< stop_date) - model_dates = model_dates[spinup_idx:stop_idx] - data = data[spinup_idx:stop_idx] - - hour_of_day = Hour.(model_dates) - mean_by_hour = [mean(data[hour_of_day .== Hour(i)]) for i in 0:23] - return mean_by_hour - end; - - # ## Experiment Setup - # - # We obtain observations from the FLUXNET site. The dataset contains multiple - # variables, but we will just use latent heat flux in this calibration. - dataset = FluxnetSimulations.get_comparison_data(site_ID, time_offset) - observations = get_diurnal_average( - (dataset.UTC_datetime, dataset.lhf), - start_date, - start_date + Day(20), - stop_date, - ); - - # Define observation noise covariance for the ensemble Kalman process. A flat - # covariance matrix is used here for simplicity. - noise_covariance = 0.05 * EKP.I; - - # ## Prior Distribution and Calibration Configuration - # - # Set up the prior distribution for the parameter and configure the ensemble - # Kalman inversion: - - # Constrained Gaussian prior for Vcmax25 with bounds [0, 1e-3], and g1 with - # bounds [0, 1000]. The first two arguments of each prior are the mean and standard - # deviation of the Gaussian function. - - priors = [ - PD.constrained_gaussian("Vcmax25", 1e-4, 5e-5, 0, 1e-3), - PD.constrained_gaussian("g1", 150, 90, 0, 1000), - ] - prior = PD.combine_distributions(priors); - - # Set the ensemble size and number of iterations - ensemble_size = 10 - N_iterations = 4; - - # ## Ensemble Kalman Inversion - # - # Initialize and run the ensemble Kalman process: - - # Sample the initial parameter ensemble from the prior distribution - initial_ensemble = EKP.construct_initial_ensemble(rng, prior, ensemble_size) - - ensemble_kalman_process = EKP.EnsembleKalmanProcess( - initial_ensemble, - observations, - noise_covariance, - EKP.Inversion(); - scheduler = EKP.DataMisfitController( - terminate_at = Inf, - on_terminate = "continue", - ), - rng, - ); - - # Run the ensemble of forward models to iteratively update the parameter - # ensemble - Logging.with_logger(SimpleLogger(devnull, Logging.Error)) do - for i in 1:N_iterations - println("Iteration $i") - params_i = EKP.get_ϕ_final(prior, ensemble_kalman_process) - G_ens = hcat([G(params_i[:, j]...) for j in 1:ensemble_size]...) - EKP.update_ensemble!(ensemble_kalman_process, G_ens) - end - end; - - # ## Results Analysis and Visualization - # Get the mean of the final parameter ensemble: - final_mean_params = EKP.get_ϕ_mean_final(prior, ensemble_kalman_process); - - ϕ = EKP.get_ϕ_final(prior, ensemble_kalman_process); - y = EKP.get_obs(ensemble_kalman_process); - metrics = EKP.get_error_metrics(ensemble_kalman_process); - - # Now, let's analyze the calibration results by examining parameter evolution - # and comparing model outputs across iterations. - # - savedir = - joinpath(pkgdir(ClimaLand), "experiments/calibration/out/$(site_ID)/out") - mkpath(savedir) - - # Save the final parameter ensemble and error metrics to CSV files for further analysis - # or record-keeping. - writedlm( - joinpath(savedir, "final_mean_parameters_$(site_ID).csv"), - final_mean_params, - ',', - ) - - writedlm( - joinpath(savedir, "final_parameter_ensemble_$(site_ID).csv"), - ϕ, - ',', - ) - - writedlm( - joinpath(savedir, "latest_data_vector_$(site_ID).csv"), - y, - ',', - ) - - df = DataFrame(metrics) - CSV.write(savedir, "error_metrics_$(site_ID).csv", df) - - # Plot the parameter ensemble evolution over iterations to visualize - # convergence: - - dim_size = sum(length.(EKP.batch(prior))) - fig = CairoMakie.Figure(size = ((dim_size + 1) * 500, 500)) - - for i in 1:dim_size - EKP.Visualize.plot_ϕ_over_iters( - fig[1, i], - ensemble_kalman_process, - prior, - i, - ) - end - - # TODO save numerical error over iterations - EKP.Visualize.plot_error_over_iters( - fig[1, dim_size + 1], - ensemble_kalman_process, - ) - CairoMakie.save("$(savedir)/constrained_params_and_error.png", fig) - # fig - - # Compare the model output between the first and last iterations to assess - # improvement: - fig = CairoMakie.Figure(size = (900, 400)) - - first_G_ensemble = EKP.get_g(ensemble_kalman_process, 1) - last_iter = EKP.get_N_iterations(ensemble_kalman_process) - last_G_ensemble = EKP.get_g(ensemble_kalman_process, last_iter) - n_ens = EKP.get_N_ens(ensemble_kalman_process) - - writedlm( - joinpath(savedir, "first_iteration_G_ensemble_$(site_ID).csv"), - first_G_ensemble, - ',', - ) - - writedlm( - joinpath(savedir, "last_iteration_G_ensemble_$(site_ID).csv"), - last_G_ensemble, - ',', - ) - - ax = Axis( - fig[1, 1]; - title = "G ensemble: first vs last iteration (n = $(n_ens), iters 1 vs $(last_iter))", - xlabel = "Observation index", - ylabel = "G", - ) - - # Plot model output of first vs last iteration ensemble - for g in eachcol(first_G_ensemble) - lines!(ax, 1:length(g), g; color = (:red, 0.6), linewidth = 1.5) - end - - for g in eachcol(last_G_ensemble) - lines!(ax, 1:length(g), g; color = (:blue, 0.6), linewidth = 1.5) - end - - # --- Add shaded noise around observations --- - xvals = 1:length(observations) - - # Say you have obs variance as `obs_var` (same length as `observations`) - # If you don’t, but only have noise covariance = 0.05*I, you can just use sqrt(0.05) - obs_std = sqrt.(0.05) .* length(observations) # or your own per-point std - - band!(ax, xvals, - observations .- 2 .* obs_std, - observations .+ 2 .* obs_std; - color = (:gray, 0.3), - transparency = true) - - # now lines on top - lines!(ax, xvals, observations; color = (:black, 0.8), linewidth = 3) - - axislegend( - ax, - [ - LineElement(color = :red, linewidth = 2), - LineElement(color = :blue, linewidth = 2), - LineElement(color = :black, linewidth = 4), - ], - ["First ensemble", "Last ensemble", "Observations"]; - position = :rb, - framevisible = false, - ) - - CairoMakie.resize_to_layout!(fig) - CairoMakie.save("$(savedir)/G_first_and_last.png", fig) - # fig -end \ No newline at end of file diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl index 3abcde3516..643b15a020 100644 --- a/experiments/integrated/fluxnet/run_fluxnet.jl +++ b/experiments/integrated/fluxnet/run_fluxnet.jl @@ -39,9 +39,9 @@ if (site_ID ∈ ("US-MOz", "US-Var", "US-NR1", "US-Ha1")) # Get the default values for this site's domain, location, and parameters (; dz_tuple, nelements, zmin, zmax) = - FluxnetSimulations.get_domain_info(FT, Val(site_ID_val)) + FluxnetSimulations.get_domain_info(FT, Val(site_ID_val)) (; time_offset, lat, long, atmos_h) = - FluxnetSimulations.get_location(FT, Val(site_ID_val)) + FluxnetSimulations.get_location(FT, Val(site_ID_val)) (; soil_ν, @@ -90,10 +90,9 @@ if (site_ID ∈ ("US-MOz", "US-Var", "US-NR1", "US-Ha1")) z0_b, ) = FluxnetSimulations.get_parameters(FT, Val(site_ID_val)) else - (; dz_tuple, nelements, zmin, zmax) = - FluxnetSimulations.get_domain_info(FT) + (; dz_tuple, nelements, zmin, zmax) = FluxnetSimulations.get_domain_info(FT) (; time_offset, lat, long, atmos_h) = - FluxnetSimulations.get_location(site_ID) + FluxnetSimulations.get_location(site_ID) end @@ -183,12 +182,8 @@ soil_domain = land_domain # ) forcing = (; atmos, radiation) -retention_parameters = (; - ν = soil_ν, - θ_r, - K_sat = soil_K_sat, - hydrology_cm = soil_hydrology_cm, -) +retention_parameters = + (; ν = soil_ν, θ_r, K_sat = soil_K_sat, hydrology_cm = soil_hydrology_cm) composition_parameters = (; ν_ss_om, ν_ss_quartz, ν_ss_gravel) soil = Soil.EnergyHydrology{FT}( diff --git a/ext/fluxnet_simulations/generic_site.jl b/ext/fluxnet_simulations/generic_site.jl index fa61bfab98..3ed6edcd3d 100644 --- a/ext/fluxnet_simulations/generic_site.jl +++ b/ext/fluxnet_simulations/generic_site.jl @@ -63,18 +63,28 @@ function FluxnetSimulations.get_parameters( soil_ν = ClimaCore.Fields.field2array(soil_params_gupta[:ν])[1], soil_K_sat = ClimaCore.Fields.field2array(soil_params_gupta[:K_sat])[1], soil_S_s = FT(1e-2), - soil_hydrology_cm = ClimaCore.MatrixFields.column_field2array(soil_params_gupta[:hydrology_cm])[1], + soil_hydrology_cm = ClimaCore.MatrixFields.column_field2array( + soil_params_gupta[:hydrology_cm], + )[1], θ_r = ClimaCore.Fields.field2array(soil_params_gupta[:θ_r])[1], - ν_ss_quartz = ClimaCore.Fields.field2array(soil_params_grids[:ν_ss_quartz])[1], + ν_ss_quartz = ClimaCore.Fields.field2array( + soil_params_grids[:ν_ss_quartz], + )[1], ν_ss_om = ClimaCore.Fields.field2array(soil_params_grids[:ν_ss_om])[1], - ν_ss_gravel = ClimaCore.Fields.field2array(soil_params_grids[:ν_ss_gravel])[1], + ν_ss_gravel = ClimaCore.Fields.field2array( + soil_params_grids[:ν_ss_gravel], + )[1], z_0m_soil = FT(0.01), z_0b_soil = FT(0.01), soil_ϵ = FT(0.98), - soil_albedo_params = ClimaLand.Soil.clm_soil_albedo_parameters(domain.space.surface), + soil_albedo_params = ClimaLand.Soil.clm_soil_albedo_parameters( + domain.space.surface, + ), soil_albedo = ClimaLand.Soil.CLMTwoBandSoilAlbedo{FT}(; - NamedTuple{keys(soil_albedo_params)}( - (ClimaCore.Fields.field2array(v)[1] for v in values(soil_albedo_params)))..., + NamedTuple{keys(soil_albedo_params)}(( + ClimaCore.Fields.field2array(v)[1] for + v in values(soil_albedo_params) + ))..., ), Ω = pft.parameters.Ω, χl = pft.parameters.χl, diff --git a/ext/fluxnet_simulations/get_fluxnet_metadata.jl b/ext/fluxnet_simulations/get_fluxnet_metadata.jl index 84a2a0bb53..c92b102533 100644 --- a/ext/fluxnet_simulations/get_fluxnet_metadata.jl +++ b/ext/fluxnet_simulations/get_fluxnet_metadata.jl @@ -63,9 +63,9 @@ function get_site_info(site_ID; fluxnet2015_metadata_path = nothing) return (; lat = nan_if_empty(site_metadata[column_name_map["latitude"]]), long = nan_if_empty(site_metadata[column_name_map["longitude"]]), - time_offset = Int64(-1 * nan_if_empty( - site_metadata[column_name_map["utc_offset"]], - )), + time_offset = Int64( + -1 * nan_if_empty(site_metadata[column_name_map["utc_offset"]]), + ), atmospheric_sensor_height = parse_array_field( site_metadata[column_name_map["atmospheric_sensor_heights"]], ), diff --git a/ext/land_sim_vis/plotting_utils.jl b/ext/land_sim_vis/plotting_utils.jl index dc817b5742..d39cf5a768 100644 --- a/ext/land_sim_vis/plotting_utils.jl +++ b/ext/land_sim_vis/plotting_utils.jl @@ -405,7 +405,7 @@ function LandSimVis.make_diurnal_timeseries( hour_of_day, data_diurnal_cycle, label = "Data", - color = "yellow", + color = "orange", ) ax.title = "$(sn): RMSD = " * diff --git a/src/Artifacts.jl b/src/Artifacts.jl index 6109661f23..073a568571 100644 --- a/src/Artifacts.jl +++ b/src/Artifacts.jl @@ -272,11 +272,14 @@ Citation: Siyan Ma, Liukang Xu, Joseph Verfaillie, Dennis Baldocchi (2023), Amer AmeriFlux CC-BY-4.0 License """ function experiment_fluxnet_data_path(site_ID; context = nothing) - try + try full_fluxnet_path = @clima_artifact("fluxnet2015", context) - dirs = filter(d -> isdir(joinpath(full_fluxnet_path, d)), readdir(full_fluxnet_path)) + dirs = filter( + d -> isdir(joinpath(full_fluxnet_path, d)), + readdir(full_fluxnet_path), + ) - match = findfirst(d -> occursin(site_ID, d), dirs) + match = findfirst(d -> occursin(site_ID, d), dirs) @assert match !== nothing "No Fluxnet data found for site ID: $site_ID" site_dir = dirs[match] @@ -289,21 +292,23 @@ function experiment_fluxnet_data_path(site_ID; context = nothing) elseif isfile(joinpath(full_fluxnet_path, site_dir, site_path_hr)) data_path = joinpath(full_fluxnet_path, site_dir, site_path_hr) else - error("There exists a directory $site_dir for site ID $site_ID, but found no data files for half-hourly or hourly data.") + error( + "There exists a directory $site_dir for site ID $site_ID, but found no data files for half-hourly or hourly data.", + ) end return data_path catch @info "Either the full fluxnet2015 dataset does not exist locally, or the site ID was not found. Falling back to the fluxnet_sites artifact which contains only US-MOz, US-Var, US-NR1, and US-Ha1." - + @assert site_ID ∈ ("US-MOz", "US-Var", "US-NR1", "US-Ha1") folder_path = @clima_artifact("fluxnet_sites", context) data_path = joinpath(folder_path, "$(site_ID).csv") return data_path end - + end function fluxnet2015_data_path(; context = nothing) diff --git a/test/integrated/fluxnet_sim.jl b/test/integrated/fluxnet_sim.jl index 866702debd..f51ef6851a 100644 --- a/test/integrated/fluxnet_sim.jl +++ b/test/integrated/fluxnet_sim.jl @@ -45,7 +45,7 @@ end @test lat == FT(42.5378) @test long == FT(-72.1715) @test atmos_h == FT(30) - + # parameters (; soil_ν, @@ -322,7 +322,7 @@ end @test nelements == 14 @test zmin == FT(-0.5) @test zmax == FT(0) - + (; time_offset, lat, long, atmos_h) = FluxnetSimulations.get_location(FT, Val(site_ID)) @@ -468,4 +468,4 @@ end @test h_leaf == FT(0.5) @test Vcmax25 == FT(2.4e-5) end -nothing \ No newline at end of file +nothing