Add plots for variables in pressure coordinates

The plots added are bias plots at 850 hPa, 500 hPa, and 250 hPa and
lat - pressure plots.

NEWS.md

@@ -21,6 +21,12 @@ by identifying where elevation is greater than 0. Note, this can lead to
 misidentification of ocean in some areas of the globe that are inland but below
 sea level (Dead Sea, Death Valley, ...).
 
+### Leaderboard for variables over longitude, latitude, time, and pressure - PR [#1094](https://github.com/CliMA/ClimaCoupler.jl/pull/1094)
+
+As a part of the post processing pipeline, bias plots for variables at the
+pressure levels of 850.0, 500.0, 250.0 hPa and bias plots over latitude and
+pressure levels are being created.
+
 
 
 ### Code cleanup

docs/src/leaderboard.md

@@ -51,3 +51,14 @@ must be initialized for each variable of interest. The CliMA model is added with
 the `RMSEVariable`. It is assumed that the `RMSEVariable` contains only the columns "DJF",
 "MAM", "JJA", "SON", and "ANN" in that order. The file `leaderboard.jl` will load the
 appropriate data into the `RMSEVariable`.
+
+### Add a new variable to compare against observations in pressure coordinates
+To add a new variable, you only need to modify the variable `sim_var_pfull_dict` in the
+function `get_sim_var_in_pfull_dict`, the variable `obs_var_dict` in the function
+`get_obs_var_in_pfull_dict`, and the variable `compare_vars_biases_plot_extrema` in the
+function `get_compare_vars_biases_plot_extrema_pfull`. The variables and functions are
+defined exactly the same as their analogous versions in the section above.
+
+It is expected that the dimensions of the variables are time, latitude, longitude, and
+pressure in no particular order and the units for the pressure dimension is expected to be
+`hPa`.

experiments/ClimaEarth/leaderboard/data_sources.jl

@@ -1,5 +1,6 @@
 import ClimaAnalysis
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
+import OrderedCollections: OrderedDict
 
 # Tuple of short names for loading simulation and observational data
 sim_obs_short_names_no_pr = [
@@ -36,7 +37,8 @@ end
 """
     get_compare_vars_biases_plot_extrema()
 
-Return a dictionary mapping short names to ranges for the bias plots.
+Return a dictionary mapping short names to ranges for the bias plots. This is
+used by the function `compute_leaderboard`.
 
 To add a variable to the leaderboard, add a key-value pair to the dictionary
 `compare_vars_biases_plot_extrema` whose key is a short name key is the same
@@ -66,7 +68,7 @@ end
     get_sim_var_dict(diagnostics_folder_path)
 
 Return a dictionary mapping short names to `OutputVar` containing preprocessed
-simulation data.
+simulation data. This is used by the function `compute_leaderboard`.
 
 To add a variable for the leaderboard, add a key-value pair to the dictionary
 `sim_var_dict` whose key is the short name of the variable and the value is an
@@ -108,7 +110,7 @@ end
     get_obs_var_dict()
 
 Return a dictionary mapping short names to `OutputVar` containing preprocessed
-observational data.
+observational data. This is used by the function `compute_leaderboard`.
 
 To add a variable for the leaderboard, add a key-value pair to the dictionary
 `obs_var_dict` whose key is the short name of the variable and the value is an
@@ -159,11 +161,11 @@ end
     get_rmse_var_dict()
 
 Return a dictionary mapping short names to `RMSEVariable` containing information about
-RMSEs of models for the short name of the variable.
+RMSEs of models for the short name of the variable. This is used by the function
+`compute_leaderboard`.
 """
 function get_rmse_var_dict()
     # Load data into RMSEVariables
-    rmse_var_names = ["pr", "rsut", "rlut"]
     rmse_var_pr = ClimaAnalysis.read_rmses(
         joinpath(@clima_artifact("cmip_model_rmse"), "pr_rmse_amip_pr_amip_5yr.csv"),
         "pr",
@@ -191,6 +193,151 @@ function get_rmse_var_dict()
     ClimaAnalysis.add_unit!(rmse_var_rlut, "CliMA", "W m^-2")
 
     # Store the rmse vars in a dictionary
-    rmse_var_dict = Dict("pr" => rmse_var_pr, "rsut" => rmse_var_rsut, "rlut" => rmse_var_rlut)
+    rmse_var_dict = OrderedDict("pr" => rmse_var_pr, "rsut" => rmse_var_rsut, "rlut" => rmse_var_rlut)
     return rmse_var_dict
 end
+
+"""
+    get_sim_var_in_pfull_dict(diagnostics_folder_path)
+
+Return a dictionary mapping short names to `OutputVar` containing preprocessed
+simulation data in pressure space. This is used by `compute_pfull_leaderboard`.
+
+To add a variable for the leaderboard, add a key-value pair to the dictionary
+`sim_var_dict` whose key is the short name of the variable and the value is an
+anonymous function that returns a `OutputVar`. For each variable, any
+preprocessing should be done in the corresponding anonymous function which
+includes unit conversion and shifting the dates.
+
+The variable should have only four dimensions: latitude, longitude, time, and
+pressure. The units of pressure should be in hPa.
+"""
+function get_sim_var_in_pfull_dict(diagnostics_folder_path)
+    simdir = ClimaAnalysis.SimDir(diagnostics_folder_path)
+    sim_var_pfull_dict = Dict{String, Any}()
+
+    short_names = ["ta", "hur", "hus"]
+    for short_name in short_names
+        sim_var_pfull_dict[short_name] =
+            () -> begin
+                sim_var = get(simdir, short_name = short_name)
+                pfull_var = get(simdir, short_name = "pfull")
+
+                (ClimaAnalysis.units(sim_var) == "") && (sim_var = ClimaAnalysis.set_units(sim_var, "unitless"))
+                (ClimaAnalysis.units(sim_var) == "kg kg^-1") &&
+                    (sim_var = ClimaAnalysis.set_units(sim_var, "unitless"))
+
+                sim_in_pfull_var = ClimaAnalysis.Atmos.to_pressure_coordinates(sim_var, pfull_var)
+                sim_in_pfull_var = ClimaAnalysis.shift_to_start_of_previous_month(sim_in_pfull_var)
+                sim_in_pfull_var = ClimaAnalysis.convert_dim_units(
+                    sim_in_pfull_var,
+                    "pfull",
+                    "hPa";
+                    conversion_function = x -> 0.01 * x,
+                )
+                return sim_in_pfull_var
+            end
+    end
+    return sim_var_pfull_dict
+end
+
+"""
+    get_obs_var_in_pfull_dict()
+
+Return a dictionary mapping short names to `OutputVar` containing preprocessed
+observational data in pressure coordinates. This is used by
+`compute_pfull_leaderboard`.
+
+To add a variable for the leaderboard, add a key-value pair to the dictionary
+`obs_var_dict` whose key is the short name of the variable and the value is an anonymous
+function that returns a `OutputVar`. The function must take in a start date
+which is used to align the times in the observational data to match the
+simulation data. The short name must be the same as in `sim_var_pfull_dict` in
+the function `get_sim_var_in_pfull_dict`. Any preprocessing is done in the
+function which includes unit conversion and shifting the dates.
+
+The variable should have only four dimensions: latitude, longitude, time, and
+pressure. The units of pressure should be in hPa.
+"""
+function get_obs_var_in_pfull_dict()
+    artifact_path = joinpath(
+        @clima_artifact("era5_monthly_averages_pressure_levels_1979_2024"),
+        "era5_monthly_averages_pressure_levels_197901-202410.nc",
+    )
+    short_names_pairs = [("ta", "t"), ("hus", "q")]
+    obs_var_dict = Dict{String, Any}()
+    for (short_name, era5_short_name) in short_names_pairs
+        obs_var_dict[short_name] =
+            (start_date) -> begin
+                obs_var = ClimaAnalysis.OutputVar(
+                    artifact_path,
+                    era5_short_name,
+                    new_start_date = start_date,
+                    shift_by = Dates.firstdayofmonth,
+                )
+                (ClimaAnalysis.units(obs_var) == "kg kg**-1") &&
+                    (obs_var = ClimaAnalysis.set_units(obs_var, "unitless"))
+                obs_var = ClimaAnalysis.Var.convert_dim_units(
+                    obs_var,
+                    "pressure_level",
+                    "hPa";
+                    conversion_function = x -> 0.01 * x,
+                )
+                return obs_var
+            end
+    end
+
+    obs_var_dict["hur"] =
+        (start_date) -> begin
+            obs_var = ClimaAnalysis.OutputVar(
+                artifact_path,
+                "r",
+                new_start_date = start_date,
+                shift_by = Dates.firstdayofmonth,
+            )
+            obs_var = ClimaAnalysis.Var.convert_dim_units(
+                obs_var,
+                "pressure_level",
+                "hPa";
+                conversion_function = x -> 0.01 * x,
+            )
+            # Convert from percentages (e.g. 120%) to decimal (1.20)
+            obs_var = ClimaAnalysis.Var.convert_units(obs_var, "unitless", conversion_function = x -> 0.01 * x)
+            return obs_var
+        end
+    return obs_var_dict
+end
+
+"""
+    get_compare_vars_biases_plot_extrema_pfull()
+
+Return a dictionary mapping short names to ranges for the bias plots. This is
+used by the function `compute_pfull_leaderboard`.
+
+To add a variable to the leaderboard, add a key-value pair to the dictionary
+`compare_vars_biases_plot_extrema` whose key is a short name key is the same
+short name in `sim_var_pfull_dict` in the function `get_sim_var_pfull_dict` and
+the value is a tuple, where the first element is the lower bound and the last
+element is the upper bound for the bias plots.
+"""
+function get_compare_vars_biases_plot_extrema_pfull()
+    compare_vars_biases_plot_extrema = Dict("ta" => (-5.0, 5.0), "hur" => (-0.4, 0.4), "hus" => (-0.0015, 0.0015))
+    return compare_vars_biases_plot_extrema
+end
+
+"""
+    get_compare_vars_biases_heatmap_extrema_pfull()
+
+Return a dictionary mapping short names to ranges for the heatmaps. This is
+used by the function `compute_pfull_leaderboard`.
+
+To add a variable to the leaderboard, add a key-value pair to the dictionary
+`compare_vars_biases_heatmap_extrema` whose key is a short name key is the same
+short name in `sim_var_pfull_dict` in the function `get_sim_var_pfull_dict` and
+the value is a tuple, where the first element is the lower bound and the last
+element is the upper bound for the bias plots.
+"""
+function get_compare_vars_biases_heatmap_extrema_pfull()
+    compare_vars_biases_heatmap_extrema = Dict("ta" => (-10.0, 10.0), "hur" => (-0.4, 0.4), "hus" => (-0.001, 0.001))
+    return compare_vars_biases_heatmap_extrema
+end