Fixed GRACE processing.

Added correct cumulative mass balance and mass balance rates, along with their
uncertainties.

analysis/analyze_scalar.py

@@ -499,6 +499,10 @@ def plot_obs_sims(
     import pism_ragis.processing  # pylint: disable=import-outside-toplevel,reimported
 
     Path(fig_dir).mkdir(exist_ok=True)
+    pdf_dir = fig_dir / Path("pdfs")
+    pdf_dir.mkdir(parents=True, exist_ok=True)
+    png_dir = fig_dir / Path("pngs")
+    png_dir.mkdir(parents=True, exist_ok=True)
 
     percentile_range = (percentiles[1] - percentiles[0]) * 100
 
@@ -648,13 +652,13 @@ def plot_obs_sims(
     prior_posterior_str = prior_str + posterior_str
 
     fig.savefig(
-        fig_dir
+        pdf_dir
         / Path(
             f"{basin}_mass_accounting_{prior_posterior_str}_filtered_by_{filtering_var}.pdf"
         )
     )
     fig.savefig(
-        fig_dir
+        png_dir
         / Path(
             f"{basin}_mass_accounting_{prior_posterior_str}_filtered_by_{filtering_var}.png",
             dpi=300,
@@ -706,6 +710,10 @@ def plot_obs_sims_3(
     import pism_ragis.processing  # pylint: disable=import-outside-toplevel,reimported
 
     Path(fig_dir).mkdir(exist_ok=True)
+    pdf_dir = fig_dir / Path("pdfs")
+    pdf_dir.mkdir(parents=True, exist_ok=True)
+    png_dir = fig_dir / Path("pngs")
+    png_dir.mkdir(parents=True, exist_ok=True)
 
     percentile_range = (percentiles[1] - percentiles[0]) * 100
 
@@ -866,70 +874,14 @@ def plot_obs_sims_3(
         posterior_str = ""
     prior_posterior_str = prior_str + posterior_str
     fig.savefig(
-        fig_dir
+        pdf_dir
         / Path(
             f"{basin}_{prior_posterior_str}_mass_accounting_filtered_by_{filtering_var}.pdf"
         )
     )
     plt.close()
 
 
-# def plot_filtered_hist():
-#     outliers_filtered_df = pd.concat([outliers_df, filtered_df]).reset_index(drop=True)
-#     # Apply the conversion function to each column
-#     outliers_filtered_df = outliers_filtered_df.apply(prp.convert_column_to_numeric)
-#     for col in ["surface.given.file", "ocean.th.file", "calving.rate_scaling.file"]:
-#         outliers_filtered_df[col] = outliers_filtered_df[col].apply(prp.simplify_path)
-#     outliers_filtered_df["surface.given.file"] = outliers_filtered_df[
-#         "surface.given.file"
-#     ].apply(prp.simplify_climate)
-#     outliers_filtered_df["ocean.th.file"] = outliers_filtered_df["ocean.th.file"].apply(
-#         prp.simplify_ocean
-#     )
-#     outliers_filtered_df["calving.rate_scaling.file"] = outliers_filtered_df[
-#         "calving.rate_scaling.file"
-#     ].apply(prp.simplify_calving)
-
-#     df = outliers_filtered_df.rename(columns=params_short_dict)
-
-#     n_params = len(params_short_dict)
-#     plt.rcParams["font.size"] = 4
-#     fig, axs = plt.subplots(
-#         5,
-#         3,
-#         sharey=True,
-#         figsize=[6.2, 6.2],
-#     )
-#     fig.subplots_adjust(hspace=1.0, wspace=0.1)
-#     for k, v in enumerate(params_short_dict.values()):
-#         legend = bool(k == 0)
-#         try:
-#             sns.histplot(
-#                 data=df,
-#                 x=v,
-#                 hue="Ensemble",
-#                 hue_order=["Filtered", "Outliers"],
-#                 palette=sim_cmap,
-#                 common_norm=False,
-#                 stat="probability",
-#                 multiple="dodge",
-#                 alpha=0.8,
-#                 linewidth=0.2,
-#                 ax=axs.ravel()[k],
-#                 legend=legend,
-#             )
-#         except:
-#             pass
-#     for ax in axs.flatten():
-#         ax.set_ylabel("")
-#         ticklabels = ax.get_xticklabels()
-#         for tick in ticklabels:
-#             tick.set_rotation(45)
-#     fn = result_dir / Path("figures") / Path("outliers_hist.pdf")
-#     fig.savefig(fn)
-#     plt.close()
-
-
 if __name__ == "__main__":
     __spec__ = None
 
@@ -1086,6 +1038,10 @@ def plot_obs_sims_3(
     result_dir = Path(options.result_dir)
     fig_dir = result_dir / Path("figures")
     fig_dir.mkdir(parents=True, exist_ok=True)
+    pdf_dir = fig_dir / Path("pdfs")
+    pdf_dir.mkdir(parents=True, exist_ok=True)
+    png_dir = fig_dir / Path("pngs")
+    png_dir.mkdir(parents=True, exist_ok=True)
 
     plt.rcParams["font.size"] = 6
 
@@ -1106,12 +1062,6 @@ def plot_obs_sims_3(
         engine=engine,
     )
 
-    # fig, ax = plt.subplots(1, 1)
-    # ds.sel(time=slice(str(filter_start_year), str(filter_end_year))).sel(
-    #     basin="GIS", ensemble_id=ensemble
-    # ).grounding_line_flux.plot(hue="exp_id", add_legend=False, ax=ax, lw=0.5)
-    # fig.savefig("grounding_line_flux_unfiltered.pdf")
-
     # Select the relevant pism_config_axis
     retreat = simulated_ds.sel(
         pism_config_axis="geometry.front_retreat.prescribed.file"
@@ -1141,7 +1091,7 @@ def plot_obs_sims_3(
     plot_outliers(
         filtered_ds.sel(basin="GIS")[outlier_variable],
         outliers_ds.sel(basin="GIS")[outlier_variable],
-        Path(fig_dir) / Path(f"{outlier_variable}_filtering.pdf"),
+        Path(pdf_dir) / Path(f"{outlier_variable}_filtering.pdf"),
     )
 
     prior_config = simulated_ds.sel(pism_config_axis=params).pism_config

data/03_prepare_mass_balance.py

@@ -34,8 +34,8 @@
 import toml
 import xarray as xr
 
+from pism_ragis.datetools import decimal_year_to_datetime
 from pism_ragis.download import download_earthaccess, download_netcdf, save_netcdf
-from pism_ragis.processing import decimal_year_to_datetime
 
 xr.set_options(keep_attrs=True)
 
@@ -95,11 +95,11 @@
     )
     for v in basin_vars:
         ds[f"cumulative_{v}"] = (ds[v] * days_in_interval).cumsum(dim="time")
-        ds[v] = ds[v].pint.to("Gt year-1")
+        ds[v] = ds[v].pint.to("Gt year^-1")
 
     for v in basin_uncertainty_vars:
         ds[f"cumulative_{v}"] = (ds[v] * days_in_interval).cumsum(dim="time")
-        ds[v] = ds[v].pint.to("Gt year-1")
+        ds[v] = ds[v].pint.to("Gt year^-1")
 
     discharge_sign = xr.DataArray(-1).pint.quantify("1")
 
@@ -124,7 +124,7 @@
         names=[
             "year",
             "cumulative_mass_balance",
-            "cumulative_mass_balance_uncertainty",
+            "mass_balance_uncertainty",
         ],
     )
 
@@ -136,15 +136,16 @@
     df["time"] = date
 
     ds = xr.Dataset.from_dataframe(df.set_index(df["time"]))
-    ds = ds.expand_dims({"basin": ["GRACE"]}, axis=-1)
+    ds["mass_balance"] = ds["cumulative_mass_balance"].diff(dim="time")
+    ds["cumulative_mass_balance_uncertainty"] = np.sqrt(
+        (ds["mass_balance_uncertainty"] ** 2).cumsum(dim="time")
+    )
     ds["cumulative_mass_balance"].attrs.update({"units": "Gt"})
     ds["cumulative_mass_balance_uncertainty"].attrs.update({"units": "Gt"})
+    ds["mass_balance"].attrs.update({"units": "Gt year^-1"})
+    ds["mass_balance_uncertainty"].attrs.update({"units": "Gt year^-1"})
+    ds = ds.expand_dims({"basin": ["GRACE"]}, axis=-1)
     fn = "grace_greenland_mass_balance.nc"
     p_fn = p / fn
     grace_ds = ds
     save_netcdf(grace_ds, p_fn)
-
-    fn = "combined_greenland_mass_balance.nc"
-    p_fn = p / fn
-    combined_ds = xr.concat([grace_ds, mankoff_ds], dim="time").sortby("time")
-    save_netcdf(combined_ds, p_fn)

pism_ragis/datetools.py

@@ -0,0 +1,75 @@
+# Copyright (C) 2023 Andy Aschwanden
+#
+# This file is part of pism-ragis.
+#
+# PISM-RAGIS is free software; you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation; either version 3 of the License, or (at your option) any later
+# version.
+#
+# PISM-RAGIS is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License
+# along with PISM; if not, write to the Free Software
+
+# pylint: disable=too-many-positional-arguments
+
+"""
+Module for date processing
+"""
+
+import datetime
+from calendar import isleap
+
+
+def decimal_year_to_datetime(decimal_year: float) -> datetime.datetime:
+    """
+    Convert a decimal year to a datetime object.
+
+    Parameters
+    ----------
+    decimal_year : float
+        The decimal year to be converted.
+
+    Returns
+    -------
+    datetime.datetime
+        The corresponding datetime object.
+
+    Notes
+    -----
+    The function calculates the date by determining the start of the year and adding
+    the fractional part of the year as days. If the resulting date has an hour value
+    of 12 or more, it rounds up to the next day and sets the time to midnight.
+    """
+    year = int(decimal_year)
+    remainder = decimal_year - year
+    start_of_year = datetime.datetime(year, 1, 1)
+    days_in_year = (datetime.datetime(year + 1, 1, 1) - start_of_year).days
+    date = start_of_year + datetime.timedelta(days=remainder * days_in_year)
+    if date.hour >= 12:
+        date = date + datetime.timedelta(days=1)
+    return date.replace(hour=0, minute=0, second=0, microsecond=0)
+
+
+def days_in_year(year: int) -> int:
+    """
+    Calculate the number of days in a given year.
+
+    Parameters
+    ----------
+    year : int
+        The year for which to calculate the number of days.
+
+    Returns
+    -------
+    int
+        The number of days in the specified year. Returns 366 if the year is a leap year, otherwise returns 365.
+    """
+    if isleap(year):
+        return 366
+    else:
+        return 365

pism_ragis/processing.py

@@ -28,7 +28,6 @@
 import re
 import shutil
 import zipfile
-from calendar import isleap
 from collections.abc import Callable
 from pathlib import Path
 from typing import Any, Dict, Hashable, List, Mapping, Union
@@ -73,56 +72,6 @@ def unzip_file(zip_path: str, extract_to: str, overwrite: bool = False) -> None:
                 zip_ref.extract(member=file, path=extract_to)
 
 
-def decimal_year_to_datetime(decimal_year: float) -> datetime.datetime:
-    """
-    Convert a decimal year to a datetime object.
-
-    Parameters
-    ----------
-    decimal_year : float
-        The decimal year to be converted.
-
-    Returns
-    -------
-    datetime.datetime
-        The corresponding datetime object.
-
-    Notes
-    -----
-    The function calculates the date by determining the start of the year and adding
-    the fractional part of the year as days. If the resulting date has an hour value
-    of 12 or more, it rounds up to the next day and sets the time to midnight.
-    """
-    year = int(decimal_year)
-    remainder = decimal_year - year
-    start_of_year = datetime.datetime(year, 1, 1)
-    days_in_year = (datetime.datetime(year + 1, 1, 1) - start_of_year).days
-    date = start_of_year + datetime.timedelta(days=remainder * days_in_year)
-    if date.hour >= 12:
-        date = date + datetime.timedelta(days=1)
-    return date.replace(hour=0, minute=0, second=0, microsecond=0)
-
-
-def days_in_year(year: int) -> int:
-    """
-    Calculate the number of days in a given year.
-
-    Parameters
-    ----------
-    year : int
-        The year for which to calculate the number of days.
-
-    Returns
-    -------
-    int
-        The number of days in the specified year. Returns 366 if the year is a leap year, otherwise returns 365.
-    """
-    if isleap(year):
-        return 366
-    else:
-        return 365
-
-
 def calculate_area(lat: np.ndarray, lon: np.ndarray) -> np.ndarray:
     """
     Calculate the area of each grid cell given arrays of latitudes and longitudes.

tests/test_datetools.py

@@ -0,0 +1,57 @@
+# Copyright (C) 2023-24 Andy Aschwanden
+#
+# This file is part of pism-ragis.
+#
+# PISM-RAGIS is free software; you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation; either version 3 of the License, or (at your option) any later
+# version.
+#
+# PISM-RAGIS is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License
+# along with PISM; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+
+"""
+Tests for date tools module.
+"""
+
+
+from pism_ragis.datetools import days_in_year
+
+
+def test_days_in_year():
+    """
+    Test the days_in_year function.
+
+    This function tests the days_in_year function for various types of years:
+    - Common year
+    - Leap year
+    - Century year that is not a leap year
+    - Century year that is a leap year
+
+    The days_in_year function is expected to return:
+    - 365 days for a common year
+    - 366 days for a leap year
+    - 365 days for a century year that is not a leap year
+    - 366 days for a century year that is a leap year
+
+    Examples
+    --------
+    >>> test_days_in_year()
+    """
+    # Test for a common year
+    assert days_in_year(2021) == 365
+
+    # Test for a leap year
+    assert days_in_year(2020) == 366
+
+    # Test for a century year that is not a leap year
+    assert days_in_year(1900) == 365
+
+    # Test for a century year that is a leap year
+    assert days_in_year(2000) == 366