diff --git a/analysis/analyze_scalar.py b/analysis/analyze_scalar.py index 53cd544..3337d66 100644 --- a/analysis/analyze_scalar.py +++ b/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,7 +874,7 @@ 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" ) @@ -874,62 +882,6 @@ def plot_obs_sims_3( 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 diff --git a/data/03_prepare_mass_balance.py b/data/03_prepare_mass_balance.py index 82b7cdc..f8c25f8 100644 --- a/data/03_prepare_mass_balance.py +++ b/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) diff --git a/pism_ragis/datetools.py b/pism_ragis/datetools.py new file mode 100644 index 0000000..b74f00b --- /dev/null +++ b/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 diff --git a/pism_ragis/processing.py b/pism_ragis/processing.py index 9a3922a..c27ac98 100644 --- a/pism_ragis/processing.py +++ b/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. diff --git a/tests/test_datetools.py b/tests/test_datetools.py new file mode 100644 index 0000000..636c529 --- /dev/null +++ b/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 diff --git a/tests/test_processing.py b/tests/test_processing.py index fe54aa7..fd124c0 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -23,40 +23,7 @@ import numpy as np import xarray as xr -from pism_ragis.processing import calculate_area, 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 +from pism_ragis.processing import calculate_area def test_drop_nonnumeric_vars():