Working ClimaLand calibration [skip ci]

Author: ph-kev

.buildkite/Manifest-v1.11.toml

@@ -429,8 +429,8 @@ weakdeps = ["GeoMakie", "Makie"]
 
 [[deps.ClimaCalibrate]]
 deps = ["Dates", "Distributed", "Distributions", "EnsembleKalmanProcesses", "JLD2", "Logging", "Random", "TOML", "YAML"]
-git-tree-sha1 = "2408ba55a30ac0f67550ef9d87cebfc372af3cca"
-repo-rev = "main"
+git-tree-sha1 = "fe9ddc24349035098116663fd0136e960d3cbb3c"
+repo-rev = "kp/identity"
 repo-url = "https://github.com/CliMA/ClimaCalibrate.jl.git"
 uuid = "4347a170-ebd6-470c-89d3-5c705c0cacc2"
 version = "0.1.2"
@@ -1251,9 +1251,9 @@ weakdeps = ["Makie"]
 
 [[deps.GeoMakie]]
 deps = ["Colors", "CoordinateTransformations", "Downloads", "GeoFormatTypes", "GeoInterface", "GeoInterfaceMakie", "GeoJSON", "Geodesy", "GeometryBasics", "GeometryOps", "ImageIO", "LinearAlgebra", "Makie", "NaturalEarth", "Proj", "Reexport", "Statistics", "StructArrays"]
-git-tree-sha1 = "8ef9756e37ca5bc892e920551f9aa683d5b8fab3"
+git-tree-sha1 = "2db1d309af35a1ad440e75b9275e4b3b4715ed39"
 uuid = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
-version = "0.7.13"
+version = "0.7.12"
 
 [[deps.Geodesy]]
 deps = ["CoordinateTransformations", "Dates", "LinearAlgebra", "StaticArrays"]

experiments/better_calibration/generate_observations.jl

@@ -145,13 +145,20 @@ function preprocess_single_era5_var(var::OutputVar, short_name, nelements)
 end
 
 if abspath(PROGRAM_FILE) == @__FILE__
+    # covar_estimator =
+    #     ClimaCalibrate.ObservationRecipe.SeasonalDiagonalCovariance(;
+    #         model_error_scale = Float32(0.05),
+    #         regularization = Float32(25.0),
+    #         ignore_nan = true,
+    #         use_latitude_weights = true,
+    #         min_cosd_lat = 0.1
+    #     )
     covar_estimator =
-        ClimaCalibrate.ObservationRecipe.SeasonalDiagonalCovariance(
-            model_error_scale = Float32(0.05),
-            regularization = Float32(25.0),
-            ignore_nan = true,
+        ClimaCalibrate.ObservationRecipe.ScalarCovariance(;
+            scalar = 25.0,
+            use_latitude_weights = true,
+            min_cosd_lat = 0.1,
         )
-
     (; nelements, sample_date_ranges, short_names) = get_config()
     @info "The number of samples is $(length(sample_date_ranges))"
 

experiments/better_calibration/getters.jl

@@ -49,7 +49,7 @@ function get_config()
         # start of DJF and SON
         sample_date_ranges = [
             (Dates.DateTime(year, 12, 1), Dates.DateTime(year + 1, 9, 1)) for
-            year in 2000:2020
+            year in 2000:2009
         ],
         spinup = Dates.Month(3), # Dates.Year(1)
         # TODO: Not sure if sample_date_ranges is the best way of determining
@@ -75,7 +75,7 @@ Return a named tuple consisting of
 function get_calibration_config()
     rng_seed = 42
     rng = Random.MersenneTwister(rng_seed)
-    return (; minibatch_size = 3, n_iterations = 7, rng = rng) # TODO Change this!
+    return (; minibatch_size = 2, n_iterations = 6, rng = rng)
 end
 
 """
@@ -119,14 +119,14 @@ end
 """
 function get_prior()
     priors = [
-        # EKP.constrained_gaussian("α_0", 0.64, 0.05, 0.2, 0.8),
         # EKP.constrained_gaussian("Δα", 0.2, 0.1, 0.0, 1.0),
         # EKP.constrained_gaussian("k", 10, 5, 2, 25),
         # EKP.constrained_gaussian("beta_snow", 0.4, 0.2, 0.1, 0.8),
         # EKP.constrained_gaussian("x0_snow", 0.4, 0.2, 0.1, 0.8),
-        EKP.constrained_gaussian("beta_snow_cover", 1.77, 0.2, 0.1, 2.5);
-        EKP.constrained_gaussian("z0_snow_cover", 0.106, 0.03, 0.0, 0.3);
-        EKP.constrained_gaussian("Δα", 0.7, 0.2, 0.0, 1.0);
+        EKP.constrained_gaussian("beta_snow_cover", 1.77, 0.2, 0.1, 2.5),
+        EKP.constrained_gaussian("z0_snow_cover", 0.106, 0.03, 0.0, 0.3),
+        EKP.constrained_gaussian("α_0", 0.7, 0.2, 0.0, 1.0),
+        # EKP.constrained_gaussian("Δα", 0.7, 0.2, 0.0, 1.0);
         # EKP.constrained_gaussian("z0_snow", 0.106, 0.05, 0.01, 0.3),
     ]
     return EKP.combine_distributions(priors)

experiments/better_calibration/model_interface.jl

@@ -74,14 +74,9 @@ function setup_model(
     calibrate_param_dict,
 )
     # TODO: This should be better when the interface improved
-    α_0 = FT(get(calibrate_param_dict, "α_0", 0.64))
-    Δα = FT(get(calibrate_param_dict, "Δα", 0.06))
-    k = FT(get(calibrate_param_dict, "k", 2))
-    beta_snow = FT(get(calibrate_param_dict, "beta_snow", 0.4))
-    x0_snow = FT(get(calibrate_param_dict, "x0_snow", 0.2))
+    α_0 = FT(get(calibrate_param_dict, "α_0", 0.7))
     beta_snow_cover = FT(get(calibrate_param_dict, "beta_snow_cover", 0.2))
     z0_snow_cover = FT(get(calibrate_param_dict, "z0_snow_cover", 0.2))
-    Δα = FT(get(calibrate_param_dict, "Δα", 0.2))
 
     time_interpolation_method =
         LONGER_RUN ? LinearInterpolation() :
@@ -305,15 +300,15 @@ function setup_model(
     )
 
     # Snow model
-    #    α_snow = Snow.ConstantAlbedoModel(FT(0.7))
+    α_snow = Snow.ConstantAlbedoModel(FT(α_0))
     # Set β = 0 in order to regain model without density dependence
-    α_snow = Snow.ZenithAngleAlbedoModel(
-        FT(α_0),
-        FT(Δα),
-        FT(k);
-        β = FT(beta_snow),
-        x0 = FT(x0_snow),
-    )
+    # α_snow = Snow.ZenithAngleAlbedoModel(
+    #     FT(α_0),
+    #     FT(Δα),
+    #     FT(k);
+    #     β = FT(beta_snow),
+    #     x0 = FT(x0_snow),
+    # )
     horz_degree_res =
         sum(ClimaLand.Domains.average_horizontal_resolution_degrees(domain)) / 2 # mean of resolution in latitude and longitude, in degrees
     scf = Snow.WuWuSnowCoverFractionModel(
@@ -356,6 +351,12 @@ function setup_model(
     return land
 end
 
+function close_diagnostics(sim)
+    for diagnostic in sim.diagnostics
+        close(diagnostic.output_writer)
+    end
+end
+
 function ClimaCalibrate.forward_model(iteration, member)
     (; output_dir, sample_date_ranges, nelements, spinup, extend) = get_config()
     ensemble_member_path =
@@ -407,5 +408,6 @@ function ClimaCalibrate.forward_model(iteration, member)
     simulation = LandSimulation(FT, start_date, stop_date, Δt, model; outdir)
     @info "Simulation is setted up. Going to solve the simulation."
     ClimaLand.Simulations.solve!(simulation)
+    close_diagnostics(simulation)
     return nothing
 end

experiments/better_calibration/observation_map.jl

@@ -90,18 +90,19 @@ function process_member_data(
     vars = map(short_names) do short_name
         var = sim_var_dict[short_name]()
         var =
-            ClimaAnalysis.average_season_across_time(var, ignore_nan = true)
+            ClimaAnalysis.average_season_across_time(var, ignore_nan = false)
 
         ocean_mask = make_ocean_mask(nelements)
         var = ocean_mask(var)
 
         # Replace all NaNs on land with the mean value
         # This is needed to stop small NaN values from stopping the calibration
-        nanmean_land_val = mean(filter(!isnan, var.data))
+        # TODO: This is a stopgap. It might be better to leave them as NaNs. As a temporary
+        # solution, we replace with 500.0
+        nanmean_land_val = 500.0 # mean(filter(!isnan, var.data))
         var = ClimaAnalysis.replace(val -> isnan(val) ? nanmean_land_val : val, var)
         var = ocean_mask(var)
 
-
         lons = ClimaAnalysis.longitudes(var)
         var = ClimaAnalysis.window(
             var,
@@ -172,16 +173,16 @@ function ClimaCalibrate.analyze_iteration(
     # EKP.TransformUnscented
     diagnostics_folder_path = joinpath(output_path, "global_diagnostics")
     try
-        # compute_monthly_leaderboard(
-        #     output_path,
-        #     diagnostics_folder_path,
-        #     "ERA5",
-        # )
-        # compute_seasonal_leaderboard(
-        #     output_path,
-        #     diagnostics_folder_path,
-        #     "ERA5",
-        # )
+        compute_monthly_leaderboard(
+            output_path,
+            diagnostics_folder_path,
+            "ERA5",
+        )
+        compute_seasonal_leaderboard(
+            output_path,
+            diagnostics_folder_path,
+            "ERA5",
+        )
     catch e
         @error "Error in `analyze_iteration`" error = e
     end

experiments/better_calibration/observation_utils.jl

@@ -83,6 +83,6 @@ function make_ocean_mask(nelements)
         landsea_var,
         NaN,
         1.0;
-        threshold = 0.99,
+        threshold = 0.5,
     )
 end