diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite.yml index 4bb4a542..2c30c1bb 100644 --- a/.github/workflows/test_suite.yml +++ b/.github/workflows/test_suite.yml @@ -73,4 +73,5 @@ jobs: python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_01_basics.py python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_02_config.py python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_03_notebooks.py - python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_04_postproc.py + python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_04_auxiliary.py + python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_05_postproc.py diff --git a/docs/calibration_options.md b/docs/calibration_options.md index a004e578..fa6117fb 100644 --- a/docs/calibration_options.md +++ b/docs/calibration_options.md @@ -6,11 +6,10 @@ Several calibration options exist, which vary with respect to complexity and com | Calibration option | Overview | Reference | | :--- | :--- | :--- | -| ['HH2015'](HH2015_target) | Finds single set of parameters.
Varies in order: $f_{snow}$, $k_{p}$, $T_{bias}$ | [Huss and Hock (2015)](https://www.frontiersin.org/articles/10.3389/feart.2015.00054/full) | -| ['HH2015mod'](HH2015mod_target) | Finds single set of parameters.
Varies in order: $k_{p}$, $T_{bias}$ | [Rounce et al. 2020](https://www.cambridge.org/core/journals/journal-of-glaciology/article/quantifying-parameter-uncertainty-in-a-largescale-glacier-evolution-model-using-bayesian-inference-application-to-high-mountain-asia/61D8956E9A6C27CC1A5AEBFCDADC0432) | -| ['emulator'](emulator_target) | Creates emulator for ['MCMC'](MCMC_target).
Finds single set of parameters with emulator following ['HH2015mod'](HH2015mod_target) | [Rounce et al. 2023](https://www.science.org/doi/10.1126/science.abo1324) | -| ['MCMC'](MCMC_target) | Finds multiple sets of parameters using Bayesian inference with [emulator](emulator_target).
Varies $f_{snow}$, $k_{p}$, $T_{bias}$ | [Rounce et al. 2023](https://www.science.org/doi/10.1126/science.abo1324) | -| ['MCMC_fullsim'](MCMC_target) | Finds multiple sets of parameters using Bayesian inference with full model simulations.
Varies $f_{snow}$, $k_{p}$, $T_{bias}$ | [Rounce et al. 2020](https://www.cambridge.org/core/journals/journal-of-glaciology/article/quantifying-parameter-uncertainty-in-a-largescale-glacier-evolution-model-using-bayesian-inference-application-to-high-mountain-asia/61D8956E9A6C27CC1A5AEBFCDADC0432) | +| ['HH2015'](HH2015_target) | Finds single set of parameters.
Varies in order: $f_{snow}$, $k_{p}$, $T_{bias}$ | [Huss and Hock, 2015](https://www.frontiersin.org/articles/10.3389/feart.2015.00054/full) | +| ['HH2015mod'](HH2015mod_target) | Finds single set of parameters.
Varies in order: $k_{p}$, $T_{bias}$ | [Rounce et al., 2020](https://www.cambridge.org/core/journals/journal-of-glaciology/article/quantifying-parameter-uncertainty-in-a-largescale-glacier-evolution-model-using-bayesian-inference-application-to-high-mountain-asia/61D8956E9A6C27CC1A5AEBFCDADC0432) | +| ['emulator'](emulator_target) | Creates emulator for ['MCMC'](MCMC_target).
Finds single set of parameters with emulator following ['HH2015mod'](HH2015mod_target) | [Rounce et al., 2023](https://www.science.org/doi/10.1126/science.abo1324) | +| ['MCMC'](MCMC_target) | Finds many sets of parameters using Bayesian inference. Setting `calib.MCMC_params.option_use_emulator=True` in ~/PyGEM/config.yaml will run Bayesian inference using the mass balance emulator. Setting `calib.MCMC_params.option_use_emulator=False` (or when performing dynamical calibration against elevation change data)$^*$ will run Bayesian inference calibration with full model simulations.
Varies $f_{snow}$, $k_{p}$, $T_{bias}$, (optionally $\rho_{ablation}$, $\rho_{accumulation}$)$^*$ | [Rounce et al., 2020](https://www.cambridge.org/core/journals/journal-of-glaciology/article/quantifying-parameter-uncertainty-in-a-largescale-glacier-evolution-model-using-bayesian-inference-application-to-high-mountain-asia/61D8956E9A6C27CC1A5AEBFCDADC0432); [2023](https://www.science.org/doi/10.1126/science.abo1324) | | [Future options](cal_custom_target) | Stay tuned for new options coming in 2023/2024! | | The output of each calibration is a .json file that holds a dictionary of the calibration options and the subsequent model parameters. Thus, the .json file will store several calibration options. Each calibration option is a key to the dictionary. The model parameters are also stored in a dictionary (i.e., a dictionary within a dictionary) with each model parameter being a key to the dictionary that provides access to a list of values for that specific model parameter. The following shows an example of how to print a list of the precipitation factors ($k_{p}$) for the calibration option specified in the input file: diff --git a/pygem/bin/postproc/postproc_binned_monthly_mass.py b/pygem/bin/postproc/postproc_binned_monthly_mass.py deleted file mode 100644 index b24b9bd2..00000000 --- a/pygem/bin/postproc/postproc_binned_monthly_mass.py +++ /dev/null @@ -1,314 +0,0 @@ -""" -Python Glacier Evolution Model (PyGEM) - -copyright © 2024 Brandon Tober David Rounce - -Distributed under the MIT license - -derive binned monthly ice thickness and mass from PyGEM simulation -""" - -# Built-in libraries -import argparse -import collections -import glob -import multiprocessing -import os -import time - -# External libraries -import numpy as np -import xarray as xr - -# pygem imports -from pygem.setup.config import ConfigManager - -# instantiate ConfigManager -config_manager = ConfigManager() -# read the config -pygem_prms = config_manager.read_config() - - -# ----- FUNCTIONS ----- -def getparser(): - """ - Use argparse to add arguments from the command line - """ - parser = argparse.ArgumentParser(description='process monthly ice thickness for PyGEM simulation') - # add arguments - parser.add_argument( - '-simpath', - action='store', - type=str, - nargs='+', - default=None, - help='path to PyGEM binned simulation (can take multiple)', - ) - parser.add_argument( - '-binned_simdir', - action='store', - type=str, - default=None, - help='directory with binned simulations for which to process monthly thickness', - ) - parser.add_argument( - '-ncores', - action='store', - type=int, - default=1, - help='number of simultaneous processes (cores) to use', - ) - - return parser - - -def get_binned_monthly(dotb_monthly, m_annual, h_annual): - """ - funciton to calculate the monthly binned ice thickness and mass - from annual climatic mass balance and annual ice thickness products - - to determine monthlyt thickness and mass, we must account for flux divergence - this is not so straight-forward, as PyGEM accounts for ice dynamics at the - end of each model year and not on a monthly timestep. - here, monthly thickness and mass is determined assuming - the flux divergence is constant throughout the year. - - annual flux divergence is first estimated by combining the annual binned change in ice - thickness and the annual binned mass balance. then, assume flux divergence is constant - throughout the year (divide annual by 12 to get monthly flux divergence). - - monthly binned flux divergence can then be combined with - monthly binned climatic mass balance to get monthly binned change in ice thickness - - - Parameters - ---------- - dotb_monthly : float - ndarray containing the climatic mass balance for each model month computed by PyGEM - shape : [#glac, #elevbins, #months] - m_annual : float - ndarray containing the average (or median) binned ice mass computed by PyGEM - shape : [#glac, #elevbins, #years] - h_annual : float - ndarray containing the average (or median) binned ice thickness at computed by PyGEM - shape : [#glac, #elevbins, #years] - - Returns - ------- - m_monthly: float - ndarray containing the binned monthly ice mass - shape : [#glac, #elevbins, #years] - h_monthly: float - ndarray containing the binned monthly ice thickness - shape : [#glac, #elevbins, #years] - """ - ### get monthly ice thickness ### - # convert mass balance from m w.e. yr^-1 to m ice yr^-1 - dotb_monthly = dotb_monthly * (pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice']) - assert dotb_monthly.shape[2] % 12 == 0, 'Number of months is not a multiple of 12!' - - # obtain annual mass balance rate, sum monthly for each year - dotb_annual = dotb_monthly.reshape(dotb_monthly.shape[0], dotb_monthly.shape[1], -1, 12).sum( - axis=-1 - ) # climatic mass balance [m ice a^-1] - - # compute the thickness change per year - delta_h_annual = np.diff(h_annual, axis=-1) # [m ice a^-1] (nbins, nyears-1) - - # compute flux divergence for each bin - flux_div_annual = dotb_annual - delta_h_annual # [m ice a^-1] - - ### to get monthly thickness and mass we require monthly flux divergence ### - # we'll assume the flux divergence is constant througohut the year (is this a good assumption?) - # ie. take annual values and divide by 12 - use numpy repeat to repeat values across 12 months - flux_div_monthly = np.repeat(flux_div_annual / 12, 12, axis=-1) - - # get monthly binned change in thickness - delta_h_monthly = dotb_monthly - flux_div_monthly # [m ice per month] - - # get binned monthly thickness = running thickness change + initial thickness - running_delta_h_monthly = np.cumsum(delta_h_monthly, axis=-1) - h_monthly = running_delta_h_monthly + h_annual[:, :, 0][:, :, np.newaxis] - - # convert to mass per unit area - m_spec_monthly = h_monthly * pygem_prms['constants']['density_ice'] - - ### get monthly mass ### - # note, binned monthly thickness and mass is currently per unit area - # obtaining binned monthly mass requires knowledge of binned glacier area - # we do not have monthly binned area (as glacier dynamics are performed on an annual timestep in PyGEM), - # so we'll resort to using the annual binned glacier mass and thickness in order to get to binned glacier area - ######################## - # first convert m_annual to bin_voluma_annual - v_annual = m_annual / pygem_prms['constants']['density_ice'] - # now get area: use numpy divide where denominator is greater than 0 to avoid divide error - # note, indexing of [:,:,1:] so that annual area array has same shape as flux_div_annual - a_annual = np.divide( - v_annual[:, :, 1:], - h_annual[:, :, 1:], - out=np.full(h_annual[:, :, 1:].shape, np.nan), - where=h_annual[:, :, 1:] > 0, - ) - - # tile to get monthly area, assuming area is constant thoughout the year - a_monthly = np.tile(a_annual, 12) - - # combine monthly thickess and area to get mass - m_monthly = m_spec_monthly * a_monthly - - return h_monthly, m_spec_monthly, m_monthly - - -def update_xrdataset(input_ds, h_monthly, m_spec_monthly, m_monthly): - """ - update xarray dataset to add new fields - - Parameters - ---------- - xrdataset : xarray Dataset - existing xarray dataset - newdata : ndarray - new data array - description: str - describing new data field - - output_ds : xarray Dataset - empty xarray dataset that contains variables and attributes to be filled in by simulation runs - encoding : dictionary - encoding used with exporting xarray dataset to netcdf - """ - # coordinates - glac_values = input_ds.glac.values - time_values = input_ds.time.values - bin_values = input_ds.bin.values - - output_coords_dict = collections.OrderedDict() - output_coords_dict['bin_thick_monthly'] = collections.OrderedDict( - [('glac', glac_values), ('bin', bin_values), ('time', time_values)] - ) - output_coords_dict['bin_mass_spec_monthly'] = collections.OrderedDict( - [('glac', glac_values), ('bin', bin_values), ('time', time_values)] - ) - output_coords_dict['bin_mass_monthly'] = collections.OrderedDict( - [('glac', glac_values), ('bin', bin_values), ('time', time_values)] - ) - - # Attributes dictionary - output_attrs_dict = {} - output_attrs_dict['bin_thick_monthly'] = { - 'long_name': 'binned monthly ice thickness', - 'units': 'm', - 'temporal_resolution': 'monthly', - 'comment': 'monthly ice thickness binned by surface elevation (assuming constant flux divergence throughout a given year)', - } - output_attrs_dict['bin_mass_spec_monthly'] = { - 'long_name': 'binned monthly specific ice mass', - 'units': 'kg m^-2', - 'temporal_resolution': 'monthly', - 'comment': 'monthly ice mass per unit area binned by surface elevation (assuming constant flux divergence throughout a given year)', - } - output_attrs_dict['bin_mass_monthly'] = { - 'long_name': 'binned monthly ice mass', - 'units': 'kg', - 'temporal_resolution': 'monthly', - 'comment': 'monthly ice mass binned by surface elevation (assuming constant flux divergence and area throughout a given year)', - } - - # Add variables to empty dataset and merge together - count_vn = 0 - encoding = {} - for vn in output_coords_dict.keys(): - empty_holder = np.zeros([len(output_coords_dict[vn][i]) for i in list(output_coords_dict[vn].keys())]) - output_ds = xr.Dataset( - {vn: (list(output_coords_dict[vn].keys()), empty_holder)}, - coords=output_coords_dict[vn], - ) - count_vn += 1 - # Merge datasets of stats into one output - if count_vn == 1: - output_ds_all = output_ds - else: - output_ds_all = xr.merge((output_ds_all, output_ds)) - # Add attributes - for vn in output_ds_all.variables: - try: - output_ds_all[vn].attrs = output_attrs_dict[vn] - except: - pass - # Encoding (specify _FillValue, offsets, etc.) - encoding[vn] = {'_FillValue': None, 'zlib': True, 'complevel': 9} - - output_ds_all['bin_thick_monthly'].values = h_monthly - output_ds_all['bin_mass_spec_monthly'].values = m_spec_monthly - output_ds_all['bin_mass_monthly'].values = m_monthly - - return output_ds_all, encoding - - -def run(simpath): - """ - create binned monthly mass change data product - Parameters - ---------- - list_packed_vars : list - list of packed variables that enable the use of parallels - Returns - ------- - binned_ds : netcdf Dataset - updated binned netcdf containing binned monthly ice thickness and mass - """ - - if os.path.isfile(simpath): - # open dataset - binned_ds = xr.open_dataset(simpath) - - # calculate monthly thickness and mass - h_monthly, m_spec_monthly, m_monthly = get_binned_monthly( - binned_ds.bin_massbalclim_monthly.values, - binned_ds.bin_mass_annual.values, - binned_ds.bin_thick_annual.values, - ) - - # update dataset to add monthly mass change - output_ds_binned, encoding_binned = update_xrdataset(binned_ds, h_monthly, m_spec_monthly, m_monthly) - - # close input ds before write - binned_ds.close() - - # append to existing binned netcdf - output_ds_binned.to_netcdf(simpath, mode='a', encoding=encoding_binned, engine='netcdf4') - - # close datasets - output_ds_binned.close() - - return - - -def main(): - time_start = time.time() - args = getparser().parse_args() - - if args.simpath: - # filter out non-file paths - simpath = [p for p in args.simpath if os.path.isfile(p)] - - elif args.binned_simdir: - # get list of sims - simpath = glob.glob(args.binned_simdir + '*.nc') - if simpath: - # number of cores for parallel processing - if args.ncores > 1: - ncores = int(np.min([len(simpath), args.ncores])) - else: - ncores = 1 - - # Parallel processing - print('Processing with ' + str(ncores) + ' cores...') - with multiprocessing.Pool(ncores) as p: - p.map(run, simpath) - - print('Total processing time:', time.time() - time_start, 's') - - -if __name__ == '__main__': - main() diff --git a/pygem/bin/postproc/postproc_binned_subannual_thick.py b/pygem/bin/postproc/postproc_binned_subannual_thick.py new file mode 100644 index 00000000..13479ae3 --- /dev/null +++ b/pygem/bin/postproc/postproc_binned_subannual_thick.py @@ -0,0 +1,356 @@ +""" +Python Glacier Evolution Model (PyGEM) + +copyright © 2024 Brandon Tober David Rounce + +Distributed under the MIT license + +derive binned subannual ice thickness and mass from PyGEM simulation +""" + +# Built-in libraries +import argparse +import collections +import glob +import json +import multiprocessing +import os +import time +from functools import partial + +# External libraries +import numpy as np +import pandas as pd +import xarray as xr + +# pygem imports +from pygem.setup.config import ConfigManager + +# instantiate ConfigManager +config_manager = ConfigManager() +# read the config +pygem_prms = config_manager.read_config() + + +# ----- FUNCTIONS ----- +def getparser(): + """ + Use argparse to add arguments from the command line + """ + parser = argparse.ArgumentParser(description='process binned subannual ice thickness for PyGEM simulation') + # add arguments + parser.add_argument( + '-simpath', + action='store', + type=str, + nargs='+', + default=None, + help='path to PyGEM binned simulation (can take multiple)', + ) + parser.add_argument( + '-simdir', + action='store', + type=str, + default=None, + help='directory with binned simulations for which to process subannual thickness', + ) + parser.add_argument( + '-ncores', + action='store', + type=int, + default=1, + help='number of simultaneous processes (cores) to use', + ) + parser.add_argument('-v', '--debug', action='store_true', help='Flag for debugging') + + return parser + + +def get_binned_subannual( + bin_massbalclim, bin_mass_annual, bin_thick_annual, dates_subannual, dates_annual, debug=False +): + """ + funciton to calculate the subannual binned ice thickness and mass + from subannual climatic mass balance and annual mass and ice thickness products. + + to determine subannual thickness and mass, we must account for flux divergence. + this is not so straight-forward, as PyGEM accounts for ice dynamics at the + end of each model year and not on a subannual timestep. + thus, subannual thickness and mass is determined assuming + the flux divergence is constant throughout the year. + + annual flux divergence is first estimated by combining the annual binned change in ice + thickness and the annual binned mass balance. then, assume flux divergence is constant + throughout the year (divide annual by the number of steps in the binned climatic mass + balance to get subannual flux divergence). + + subannual binned flux divergence can then be combined with + subannual binned climatic mass balance to get subannual binned change in ice thickness and mass. + + + Parameters + ---------- + bin_massbalclim : ndarray + climatic mass balance [m w.e. yr^-1] with subannual timesteps (monthly/daily) + shape : [#glac, #elevbins, #steps] + bin_mass_annual : ndarray + annual binned ice mass computed by PyGEM [kg] + shape : [#glac, #elevbins, #years] + bin_thick_annual : ndarray + annual binned glacier thickness [m ice] + shape : [#glac, #elevbins, #years] + dates_subannual : array-like of datetime-like + dates associated with `bin_massbalclim` (subannual) + dates_annual : array-like of datetime-like + dates associated with `bin_thick_annual` and `bin_mass_annual` (annual, values correspond to start of the year) + + + Returns + ------- + h_subannual : ndarray + subannual binned ice thickness [m ice] + shape : [#glac, #elevbins, #steps] + m_spec_subannual : ndarray + subannual binned specific ice mass [kg m^-2] + shape : [#glac, #elevbins, #steps] + m_subannual : ndarray + subannual binned glacier mass [kg] + shape : [#glac, #elevbins, #steps] + """ + + n_glac, n_bins, n_steps = bin_massbalclim.shape + years_annual = np.array([d.year for d in dates_annual]) + years_subannual = np.array([d.year for d in dates_subannual]) + yrs = np.unique(years_subannual) + nyrs = len(yrs) + assert nyrs > 1, 'Need at least two annual steps for flux divergence estimation' + + # --- Step 1: convert mass balance from m w.e. to m ice --- + rho_w = pygem_prms['constants']['density_water'] + rho_i = pygem_prms['constants']['density_ice'] + bin_massbalclim_ice = bin_massbalclim * (rho_w / rho_i) + + # --- Step 2: compute annual cumulative mass balance --- + # Initialize annual cumulative mass balance (exclude last year for flux calculation) + bin_massbalclim_annual = np.zeros((n_glac, n_bins, nyrs)) + for i, year in enumerate(yrs): + idx = np.where(years_subannual == year)[0] + bin_massbalclim_annual[:, :, i] = bin_massbalclim_ice[:, :, idx].sum(axis=-1) + + # --- Step 3: compute annual thickness change --- + bin_delta_thick_annual = np.diff(bin_thick_annual, axis=-1) # [m ice yr^-1] + + # --- Step 4: compute annual flux divergence --- + bin_flux_divergence_annual = bin_massbalclim_annual - bin_delta_thick_annual # [m ice yr^-1] + + # --- Step 5: expand flux divergence to subannual steps --- + bin_flux_divergence_subannual = np.zeros_like(bin_massbalclim_ice) + for i, year in enumerate(yrs): + idx = np.where(years_subannual == year)[0] + bin_flux_divergence_subannual[:, :, idx] = bin_flux_divergence_annual[:, :, i][:, :, np.newaxis] / len(idx) + + # --- Step 6: compute subannual thickness change --- + bin_delta_thick_subannual = bin_massbalclim_ice - bin_flux_divergence_subannual + + # --- Step 7: calculate subannual thickness --- + running_bin_delta_thick_subannual = np.cumsum(bin_delta_thick_subannual, axis=-1) + bin_thick_subannual = running_bin_delta_thick_subannual + bin_thick_annual[:, :, 0][:, :, np.newaxis] + + # --- Step 8: compute glacier volume and area on subannual timestep --- + bin_volume_annual = bin_mass_annual / rho_i # annual volume [m^3] per bin + bin_area_annual = np.divide( + bin_volume_annual[:, :, 1:], # exclude first year to match flux_div_annual + bin_thick_annual[:, :, 1:], + out=np.full(bin_thick_annual[:, :, 1:].shape, np.nan), + where=bin_thick_annual[:, :, 1:] > 0, + ) + + # --- Step 9 : compute subannual glacier mass --- + # First expand area to subannual steps + bin_area_subannual = np.full(bin_massbalclim_ice.shape, np.nan) + for i, year in enumerate(yrs): + idx = np.where(years_subannual == year)[0] + bin_area_subannual[:, :, idx] = bin_area_annual[:, :, i][:, :, np.newaxis] + + # multiply by ice density to get subannual mass + bin_mass_subannual = bin_thick_subannual * rho_i * bin_area_subannual + + # --- Step 10: debug check --- + if debug: + for i, year in enumerate(yrs): + # get last subannual index of that year + idx = np.where(years_subannual == year)[0][-1] + diff = bin_thick_subannual[:, :, idx] - bin_thick_annual[:, :, i + 1] + print(f'Year {year}, subannual idx: {idx}') + print('Max diff:', np.max(np.abs(diff))) + print('Min diff:', np.min(np.abs(diff))) + print('Mean diff:', np.mean(diff)) + print() + # optional assertion + np.testing.assert_allclose( + bin_thick_subannual[:, :, idx], + bin_thick_annual[:, :, i + 1], + rtol=1e-6, + atol=1e-12, + err_msg=f'Mismatch in thickness for year {year}', + ) + + return bin_thick_subannual, bin_mass_subannual + + +def update_xrdataset(input_ds, bin_thick, bin_mass, timestep): + """ + update xarray dataset to add new fields + + Parameters + ---------- + xrdataset : xarray Dataset + existing xarray dataset + newdata : ndarray + new data array + description: str + describing new data field + + output_ds : xarray Dataset + empty xarray dataset that contains variables and attributes to be filled in by simulation runs + encoding : dictionary + encoding used with exporting xarray dataset to netcdf + """ + # coordinates + glac_values = input_ds.glac.values + time_values = input_ds.time.values + bin_values = input_ds.bin.values + + output_coords_dict = collections.OrderedDict() + output_coords_dict['bin_thick'] = collections.OrderedDict( + [('glac', glac_values), ('bin', bin_values), ('time', time_values)] + ) + output_coords_dict['bin_mass'] = collections.OrderedDict( + [('glac', glac_values), ('bin', bin_values), ('time', time_values)] + ) + + # Attributes dictionary + output_attrs_dict = {} + output_attrs_dict['bin_thick'] = { + 'long_name': 'binned ice thickness', + 'units': 'm', + 'temporal_resolution': timestep, + 'comment': 'subannual ice thickness binned by surface elevation (assuming constant flux divergence throughout a given year)', + } + output_attrs_dict['bin_mass'] = { + 'long_name': 'binned ice mass', + 'units': 'kg', + 'temporal_resolution': timestep, + 'comment': 'subannual ice mass binned by surface elevation (assuming constant flux divergence and area throughout a given year)', + } + + # Add variables to empty dataset and merge together + count_vn = 0 + encoding = {} + for vn in output_coords_dict.keys(): + empty_holder = np.zeros([len(output_coords_dict[vn][i]) for i in list(output_coords_dict[vn].keys())]) + output_ds = xr.Dataset( + {vn: (list(output_coords_dict[vn].keys()), empty_holder)}, + coords=output_coords_dict[vn], + ) + count_vn += 1 + # Merge datasets of stats into one output + if count_vn == 1: + output_ds_all = output_ds + else: + output_ds_all = xr.merge((output_ds_all, output_ds)) + # Add attributes + for vn in output_ds_all.variables: + try: + output_ds_all[vn].attrs = output_attrs_dict[vn] + except: + pass + # Encoding (specify _FillValue, offsets, etc.) + encoding[vn] = {'_FillValue': None, 'zlib': True, 'complevel': 9} + + output_ds_all['bin_thick'].values = bin_thick + output_ds_all['bin_mass'].values = bin_mass + + return output_ds_all, encoding + + +def run(simpath, debug): + """ + create binned subannual mass change data product + Parameters + ---------- + list_packed_vars : list + list of packed variables that enable the use of parallels + Returns + ------- + binned_ds : netcdf Dataset + updated binned netcdf containing binned subannual ice thickness and mass + """ + + if os.path.isfile(simpath): + # open dataset + binned_ds = xr.open_dataset(simpath) + + # get model time tables + timestep = json.loads(binned_ds.attrs['model_parameters'])['timestep'] + # get model dates + dates_annual = pd.to_datetime([f'{y}-01-01' for y in binned_ds.year.values]) + dates_subannual = pd.to_datetime(binned_ds.time.values) + + # calculate subannual thickness and mass + bin_thick, bin_mass = get_binned_subannual( + bin_massbalclim=binned_ds.bin_massbalclim.values, + bin_mass_annual=binned_ds.bin_mass_annual.values, + bin_thick_annual=binned_ds.bin_thick_annual.values, + dates_subannual=dates_subannual, + dates_annual=dates_annual, + debug=debug, + ) + + # update dataset to add subannual binned thickness and mass + output_ds_binned, encoding_binned = update_xrdataset( + binned_ds, bin_thick=bin_thick, bin_mass=bin_mass, timestep=timestep + ) + + # close input ds before write + binned_ds.close() + + # append to existing binned netcdf + output_ds_binned.to_netcdf(simpath, mode='a', encoding=encoding_binned, engine='netcdf4') + + # close datasets + output_ds_binned.close() + + return + + +def main(): + time_start = time.time() + args = getparser().parse_args() + + if args.simpath: + # filter out non-file paths + simpath = [p for p in args.simpath if os.path.isfile(p)] + elif args.simdir: + # get list of sims + simpath = sorted(glob.glob(args.simdir + '/*.nc')) + if simpath: + # number of cores for parallel processing + if args.ncores > 1: + ncores = int(np.min([len(simpath), args.ncores])) + else: + ncores = 1 + + # set up partial function with debug argument + run_partial = partial(run, debug=args.debug) + + # Parallel processing + print('Processing with ' + str(ncores) + ' cores...') + with multiprocessing.Pool(ncores) as p: + p.map(run_partial, simpath) + + print('Total processing time:', time.time() - time_start, 's') + + +if __name__ == '__main__': + main() diff --git a/pygem/bin/postproc/postproc_compile_simulations.py b/pygem/bin/postproc/postproc_compile_simulations.py index 391fff75..2247687e 100644 --- a/pygem/bin/postproc/postproc_compile_simulations.py +++ b/pygem/bin/postproc/postproc_compile_simulations.py @@ -11,6 +11,7 @@ # imports import argparse import glob +import json import multiprocessing import os import time @@ -57,67 +58,67 @@ # define metadata for each variable var_metadata = { - 'glac_runoff_monthly': { + 'glac_runoff': { 'long_name': 'glacier-wide runoff', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', 'comment': 'runoff from the glacier terminus, which moves over time', }, - 'offglac_runoff_monthly': { + 'offglac_runoff': { 'long_name': 'off-glacier-wide runoff', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', 'comment': 'off-glacier runoff from area where glacier no longer exists', }, - 'glac_acc_monthly': { + 'glac_acc': { 'long_name': 'glacier-wide accumulation, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', 'comment': 'only the solid precipitation', }, - 'glac_melt_monthly': { + 'glac_melt': { 'long_name': 'glacier-wide melt, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', }, - 'glac_refreeze_monthly': { + 'glac_refreeze': { 'long_name': 'glacier-wide refreeze, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', }, - 'glac_frontalablation_monthly': { + 'glac_frontalablation': { 'long_name': 'glacier-wide frontal ablation, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', 'comment': ( 'mass losses from calving, subaerial frontal melting, sublimation above the waterline and ' 'subaqueous frontal melting below the waterline; positive values indicate mass lost like melt' ), }, - 'glac_snowline_monthly': { + 'glac_snowline': { 'long_name': 'transient snowline altitude above mean sea level', 'units': 'm', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', 'comment': 'transient snowline is altitude separating snow from ice/firn', }, - 'glac_massbaltotal_monthly': { + 'glac_massbaltotal': { 'long_name': 'glacier-wide total mass balance, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', 'comment': 'total mass balance is the sum of the climatic mass balance and frontal ablation', }, - 'glac_prec_monthly': { + 'glac_prec': { 'long_name': 'glacier-wide precipitation (liquid)', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', 'comment': 'only the liquid precipitation, solid precipitation excluded', }, - 'glac_mass_monthly': { + 'glac_mass': { 'long_name': 'glacier mass', 'units': 'kg', - 'temporal_resolution': 'monthly', + 'temporal_resolution': '', 'comment': ( - 'mass of ice based on area and ice thickness at start of the year and the monthly total mass balance' + 'mass of ice based on area and ice thickness at start of the year and the total mass balance over during the model year' ), }, 'glac_area_annual': { @@ -219,13 +220,16 @@ def run(args): + '/' + sim_climate_scenario + '/stats/' - + f'*{gcm}_{sim_climate_scenario}_{realizations[0]}_{calibration}_ba{bias_adj}_*_{sim_startyear}_{sim_endyear}_all.nc'.replace( + + f'*{gcm}_{sim_climate_scenario}_{realizations[0]}_{calibration}_ba{bias_adj}_*sets_{sim_startyear}_{sim_endyear}_all.nc'.replace( '__', '_' ) )[0] else: fp = glob.glob( - base_dir + gcm + '/stats/' + f'*{gcm}_{calibration}_ba{bias_adj}_*_{sim_startyear}_{sim_endyear}_all.nc' + base_dir + + gcm + + '/stats/' + + f'*{gcm}_{calibration}_ba{bias_adj}_*sets_{sim_startyear}_{sim_endyear}_all.nc' )[0] # get number of sets from file name nsets = fp.split('/')[-1].split('_')[-4] @@ -233,6 +237,12 @@ def run(args): ds_glac = xr.open_dataset(fp) year_values = ds_glac.year.values time_values = ds_glac.time.values + # get model time step + timestep = json.loads(ds_glac.attrs['model_parameters'])['timestep'] + if timestep.lower() in ['daily']: + period_name = 'day' + elif timestep.lower() in ['monthly']: + period_name = 'month' # check if desired vars are in ds ds_vars = list(ds_glac.keys()) missing_vars = list(set(vars) - set(ds_vars)) @@ -317,19 +327,19 @@ def run(args): try: arr = getattr(ds_glac, var).values except AttributeError: - if 'monthly' in var: - arr = np.full((1, len(time_values)), np.nan) - elif 'annual' in var: + if 'annual' in var: arr = np.full((1, year_values.shape[0]), np.nan) + else: + arr = np.full((1, len(time_values)), np.nan) reg_gcm_data[var].append(arr) # if glacier output DNE in sim output file, create empty nan arrays to keep record of missing glaciers except: for var in vars: - if 'monthly' in var: - arr = np.full((1, len(time_values)), np.nan) - elif 'annual' in var: + if 'annual' in var: arr = np.full((1, year_values.shape[0]), np.nan) + else: + arr = np.full((1, len(time_values)), np.nan) reg_gcm_data[var].append(arr) # stack all individual glacier data for each var, for the current GCM and append as 2D array to list of reg_all_gcms_data[var] @@ -408,6 +418,8 @@ def run(args): ]: if attr_key in meta: ds[var].attrs[attr_key] = meta[attr_key] + if attr_key == 'temporal_resolution' and 'annual' not in var: + ds[var].attrs[attr_key] = timestep ds[var].attrs['grid_mapping'] = 'crs' # crs attributes - same for all vars @@ -425,10 +437,9 @@ def run(args): if 'annual' in var: ds.time.attrs['range'] = str(year_values[0]) + ' - ' + str(year_values[-1]) ds.time.attrs['comment'] = 'years referring to the start of each year' - elif 'monthly' in var: + else: ds.time.attrs['range'] = str(time_values[0]) + ' - ' + str(time_values[-1]) - ds.time.attrs['comment'] = 'start of the month' - + ds.time.attrs['comment'] = f'start of the {period_name}' ds.RGIId.attrs['long_name'] = 'Randolph Glacier Inventory Id' ds.RGIId.attrs['comment'] = 'RGIv6.0 (https://nsidc.org/data/nsidc-0770/versions/6)' ds.RGIId.attrs['cf_role'] = 'timeseries_id' @@ -590,18 +601,18 @@ def main(): parser.add_argument( '-vars', type=str, - help='comm delimited list of PyGEM variables to compile (can take multiple, ex. "monthly_mass annual_area")', + help='comm delimited list of PyGEM variables to compile (can take multiple, ex. "glac_mass glac_area_annual")', choices=[ - 'glac_runoff_monthly', - 'offglac_runoff_monthly', - 'glac_acc_monthly', - 'glac_melt_monthly', - 'glac_refreeze_monthly', - 'glac_frontalablation_monthly', - 'glac_snowline_monthly', - 'glac_massbaltotal_monthly', - 'glac_prec_monthly', - 'glac_mass_monthly', + 'glac_runoff', + 'offglac_runoff', + 'glac_acc', + 'glac_melt', + 'glac_refreeze', + 'glac_frontalablation', + 'glac_snowline', + 'glac_massbaltotal', + 'glac_prec', + 'glac_mass', 'glac_mass_annual', 'glac_area_annual', 'glac_ELA_annual', @@ -669,16 +680,16 @@ def main(): if not vars: vars = [ - 'glac_runoff_monthly', - 'offglac_runoff_monthly', - 'glac_acc_monthly', - 'glac_melt_monthly', - 'glac_refreeze_monthly', - 'glac_frontalablation_monthly', - 'glac_snowline_monthly', - 'glac_massbaltotal_monthly', - 'glac_prec_monthly', - 'glac_mass_monthly', + 'glac_runoff', + 'offglac_runoff', + 'glac_acc', + 'glac_melt', + 'glac_refreeze', + 'glac_frontalablation', + 'glac_snowline', + 'glac_massbaltotal', + 'glac_prec', + 'glac_mass', 'glac_mass_annual', 'glac_area_annual', 'glac_ELA_annual', diff --git a/pygem/bin/postproc/postproc_monthly_mass.py b/pygem/bin/postproc/postproc_subannual_mass.py similarity index 52% rename from pygem/bin/postproc/postproc_monthly_mass.py rename to pygem/bin/postproc/postproc_subannual_mass.py index 819c2591..7931d937 100644 --- a/pygem/bin/postproc/postproc_monthly_mass.py +++ b/pygem/bin/postproc/postproc_subannual_mass.py @@ -5,18 +5,22 @@ Distributed under the MIT license -derive monthly glacierwide mass for PyGEM simulation using annual glacier mass and monthly total mass balance +derive sub-annual glacierwide mass for PyGEM simulation using annual glacier mass and sub-annual total mass balance """ # Built-in libraries import argparse import collections import glob +import json import multiprocessing import os import time +from functools import partial +import matplotlib.pyplot as plt import numpy as np +import pandas as pd # External libraries import xarray as xr @@ -36,7 +40,7 @@ def getparser(): Use argparse to add arguments from the command line """ parser = argparse.ArgumentParser( - description='process monthly glacierwide mass from annual mass and total monthly mass balance' + description='process sub-annual glacierwide mass from annual mass and total sub-annual mass balance' ) # add arguments parser.add_argument( @@ -51,7 +55,7 @@ def getparser(): action='store', type=str, default=None, - help='directory with glacierwide simulation outputs for which to process monthly mass', + help='directory with glacierwide simulation outputs for which to process sub-annual mass', ) parser.add_argument( '-ncores', @@ -60,14 +64,14 @@ def getparser(): default=1, help='number of simultaneous processes (cores) to use', ) - + parser.add_argument('-v', '--debug', action='store_true', help='Flag for debugging') return parser -def get_monthly_mass(glac_mass_annual, glac_massbaltotal_monthly): +def get_subannual_mass(df_annual, df_sub, debug=False): """ - funciton to calculate the monthly glacier mass - from annual glacier mass and monthly total mass balance + funciton to calculate the sub-annual glacier mass + from annual glacier mass and sub-annual total mass balance Parameters ---------- @@ -75,35 +79,65 @@ def get_monthly_mass(glac_mass_annual, glac_massbaltotal_monthly): ndarray containing the annual glacier mass for each year computed by PyGEM shape: [#glac, #years] unit: kg - glac_massbaltotal_monthly : float - ndarray containing the monthly total mass balance computed by PyGEM - shape: [#glac, #months] + glac_massbaltotal : float + ndarray containing the total mass balance computed by PyGEM + shape: [#glac, #steps] unit: kg Returns ------- - glac_mass_monthly: float - ndarray containing the monthly glacier mass - shape : [#glac, #months] + glac_mass: float + ndarray containing the running glacier mass + shape : [#glac, #steps] unit: kg """ - # get running total monthly mass balance - reshape into subarrays of all values for a given year, then take cumulative sum - oshape = glac_massbaltotal_monthly.shape - running_glac_massbaltotal_monthly = ( - np.reshape(glac_massbaltotal_monthly, (-1, 12), order='C').cumsum(axis=-1).reshape(oshape) - ) - # tile annual mass to then superimpose atop running glacier mass balance (trim off final year from annual mass) - glac_mass_monthly = np.repeat(glac_mass_annual[:, :-1], 12, axis=-1) + # ensure datetime and sorted + df_annual['time'] = pd.to_datetime(df_annual['time']) + df_sub['time'] = pd.to_datetime(df_sub['time']) + df_annual = df_annual.sort_values('time').reset_index(drop=True) + df_sub = df_sub.sort_values('time').reset_index(drop=True) + + # year columns + df_annual['year'] = df_annual['time'].dt.year + df_sub['year'] = df_sub['time'].dt.year + + # map annual starting mass to sub rows + annual_by_year = df_annual.set_index('year')['mass'] + df_sub['annual_mass'] = df_sub['year'].map(annual_by_year) + + # shift massbaltotal within each year so the Jan value doesn't affect Jan mass itself + # i.e., massbaltotal at Jan-01 contributes to Feb-01 mass + df_sub['mb_shifted'] = df_sub.groupby('year')['massbaltotal'].shift(1).fillna(0.0) + + # cumulative sum of shifted values within each year + df_sub['cum_mb_since_year_start'] = df_sub.groupby('year')['mb_shifted'].cumsum() + + # compute sub-annual mass + df_sub['mass'] = df_sub['annual_mass'] + df_sub['cum_mb_since_year_start'] - # add annual mass values to running glacier mass balance - glac_mass_monthly += running_glac_massbaltotal_monthly + if debug: + # --- Quick plot of Jan start points (sub vs annual) --- + # Plot all sub-annual masses as a line + plt.figure(figsize=(12, 5)) + plt.plot(df_sub['time'], df_sub['mass'], label='Sub-annual mass', color='blue') - return glac_mass_monthly + # Overlay annual masses as points/line + plt.plot(df_annual['time'], df_annual['mass'], 'o--', label='Annual mass', color='orange', markersize=6) + # Labels and legend + plt.xlabel('Time') + plt.ylabel('Glacier Mass') + plt.title('Sub-annual Glacier Mass vs Annual Mass') + plt.legend() + plt.tight_layout() + plt.show() -def update_xrdataset(input_ds, glac_mass_monthly): + return df_sub['mass'].values + + +def update_xrdataset(input_ds, glac_mass, timestep): """ update xarray dataset to add new fields @@ -126,15 +160,15 @@ def update_xrdataset(input_ds, glac_mass_monthly): time_values = input_ds.time.values output_coords_dict = collections.OrderedDict() - output_coords_dict['glac_mass_monthly'] = collections.OrderedDict([('glac', glac_values), ('time', time_values)]) + output_coords_dict['glac_mass'] = collections.OrderedDict([('glac', glac_values), ('time', time_values)]) # Attributes dictionary output_attrs_dict = {} - output_attrs_dict['glac_mass_monthly'] = { + output_attrs_dict['glac_mass'] = { 'long_name': 'glacier mass', 'units': 'kg', - 'temporal_resolution': 'monthly', - 'comment': 'monthly glacier mass', + 'temporal_resolution': timestep, + 'comment': 'glacier mass', } # Add variables to empty dataset and merge together @@ -160,15 +194,14 @@ def update_xrdataset(input_ds, glac_mass_monthly): pass # Encoding (specify _FillValue, offsets, etc.) encoding[vn] = {'_FillValue': None, 'zlib': True, 'complevel': 9} - - output_ds_all['glac_mass_monthly'].values = glac_mass_monthly + output_ds_all['glac_mass'].values = glac_mass[np.newaxis, :] return output_ds_all, encoding -def run(simpath): +def run(simpath, debug=False): """ - create monthly mass data product + create sub-annual mass data product Parameters ---------- simpath : str @@ -178,16 +211,27 @@ def run(simpath): try: # open dataset statsds = xr.open_dataset(simpath) - - # calculate monthly mass - pygem glac_massbaltotal_monthly is in units of m3, so convert to mass using density of ice - glac_mass_monthly = get_monthly_mass( - statsds.glac_mass_annual.values, - statsds.glac_massbaltotal_monthly.values * pygem_prms['constants']['density_ice'], + timestep = json.loads(statsds.attrs['model_parameters'])['timestep'] + yvals = statsds.year.values + # convert to pandas dataframe with annual mass + annual_df = pd.DataFrame( + {'time': pd.to_datetime([f'{y}-01-01' for y in yvals]), 'mass': statsds.glac_mass_annual[0].values} ) + tvals = statsds.time.values + # convert to pandas dataframe with sub-annual mass balance + steps_df = pd.DataFrame( + { + 'time': pd.to_datetime([t.strftime('%Y-%m-%d') for t in tvals]), + 'massbaltotal': statsds.glac_massbaltotal[0].values * pygem_prms['constants']['density_ice'], + } + ) + + # calculate sub-annual mass - pygem glac_massbaltotal is in units of m3, so convert to mass using density of ice + glac_mass = get_subannual_mass(annual_df, steps_df, debug=debug) statsds.close() - # update dataset to add monthly mass change - output_ds_stats, encoding = update_xrdataset(statsds, glac_mass_monthly) + # update dataset to add sub-annual mass change + output_ds_stats, encoding = update_xrdataset(statsds, glac_mass, timestep) # close input ds before write statsds.close() @@ -225,10 +269,13 @@ def main(): else: ncores = 1 + # set up partial function with debug argument + run_partial = partial(run, debug=args.debug) + # Parallel processing print('Processing with ' + str(ncores) + ' cores...') with multiprocessing.Pool(ncores) as p: - p.map(run, simpath) + p.map(run_partial, simpath) print('Total processing time:', time.time() - time_start, 's') diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py index 180a458d..d7a3faa6 100755 --- a/pygem/bin/run/run_calibration.py +++ b/pygem/bin/run/run_calibration.py @@ -318,7 +318,7 @@ def calculate_elev_change_1d( ev_model.mb_model.glac_wide_massbaltotal + ev_model.mb_model.glac_wide_frontalablation ) - mod_glacierwide_mb_mwea = ( + glacierwide_mb_mwea = ( mbmod.glac_wide_massbaltotal[gdir.mbdata['t1_idx'] : gdir.mbdata['t2_idx'] + 1].sum() / mbmod.glac_wide_area_annual[0] / gdir.mbdata['nyears'] @@ -329,62 +329,63 @@ def calculate_elev_change_1d( except RuntimeError: return float('-inf'), float('-inf') - ### get monthly ice thickness - # grab components of interest - thickness_m = ds[0].thickness_m.values.T # glacier thickness [m ice], (nbins, nyears) - - # set any < 0 thickness to nan - thickness_m[thickness_m <= 0] = np.nan - - # climatic mass balance - dotb_monthly = mbmod.glac_bin_massbalclim # climatic mass balance [m w.e.] per month - # convert to m ice - dotb_monthly = dotb_monthly * (pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice']) - - ### to get monthly thickness and mass we require monthly flux divergence ### - # we'll assume the flux divergence is constant througohut the year - # ie. take annual values and divide by 12 - use numpy repeat to repeat values across 12 months - flux_div_monthly_mmo = np.repeat(-ds[0].flux_divergence_myr.values.T[:, 1:] / 12, 12, axis=-1) - - # get monthly binned change in thickness - delta_h_monthly = dotb_monthly - flux_div_monthly_mmo # [m ice per month] - - # get binned monthly thickness = running thickness change + initial thickness - running_delta_h_monthly = np.cumsum(delta_h_monthly, axis=-1) - h_monthly = running_delta_h_monthly + thickness_m[:, 0][:, np.newaxis] - + ### get subannual elevation change + + # --- Step 1: convert mass balance from m w.e. to m ice --- + rho_w = pygem_prms['constants']['density_water'] + rho_i = pygem_prms['constants']['density_ice'] + bin_massbalclim_ice = mbmod.glac_bin_massbalclim * (rho_w / rho_i) # binned climatic mass balance (nbins, nsteps) + + # --- Step 2: expand flux divergence to subannual steps --- + # assume the flux divergence is constant througohut the year + # ie. take annual values and divide spread uniformly throughout model year + bin_flux_divergence_subannual = np.zeros_like(bin_massbalclim_ice) + for i, year in enumerate(gdir.dates_table.year.unique()): + idx = np.where(gdir.dates_table.year.values == year)[0] + bin_flux_divergence_subannual[:, idx] = -ds[0].flux_divergence_myr.values.T[:, i + 1][:, np.newaxis] / len( + idx + ) # note, oggm flux_divergence_myr is opposite sign of convention, hence negative + + # --- Step 3: compute subannual thickness change --- + bin_delta_thick_subannual = bin_massbalclim_ice - bin_flux_divergence_subannual # [m ice] + + # --- Step 4: calculate subannual thickness --- + # calculate binned subannual thickness = running thickness change + initial thickness + bin_thick_initial = ds[0].thickness_m.isel(time=0).values # initial glacier thickness [m ice], (nbins) + running_bin_delta_thick_subannual = np.cumsum(bin_delta_thick_subannual, axis=-1) + bin_thick_subannual = running_bin_delta_thick_subannual + bin_thick_initial[:, np.newaxis] + + # --- Step 5: rebin subannual thickness --- # get surface height at the specified reference year - ref_surface_h = ds[0].bed_h.values + ds[0].thickness_m.sel(time=gdir.elev_change_1d['ref_dem_year']).values - + ref_surface_height = ds[0].bed_h.values + ds[0].thickness_m.sel(time=gdir.elev_change_1d['ref_dem_year']).values # aggregate model bin thicknesses as desired with warnings.catch_warnings(): warnings.filterwarnings('ignore') - h_monthly = np.column_stack( + bin_thick_subannual = np.column_stack( [ stats.binned_statistic( - x=ref_surface_h, + x=ref_surface_height, values=x, statistic=np.nanmean, bins=gdir.elev_change_1d['bin_edges'], )[0] - for x in h_monthly.T + for x in bin_thick_subannual.T ] ) - # interpolate over any empty bins - h_monthly_ = np.column_stack([interp1d_fill_gaps(x.copy()) for x in h_monthly.T]) + bin_thick_subannual = np.column_stack([interp1d_fill_gaps(x.copy()) for x in bin_thick_subannual.T]) - # difference each set of inds in diff_inds_map - mod_elev_change_1d = np.column_stack( + # --- Step 6: calculate elevation change --- + elev_change_1d = np.column_stack( [ - h_monthly_[:, tup[1]] - h_monthly_[:, tup[0]] + bin_thick_subannual[:, tup[1]] - bin_thick_subannual[:, tup[0]] if tup[0] is not None and tup[1] is not None - else np.full(h_monthly_.shape[0], np.nan) + else np.full(bin_thick_subannual.shape[0], np.nan) for tup in gdir.elev_change_1d['model2obs_inds_map'] ] ) - return mod_glacierwide_mb_mwea, mod_elev_change_1d + return glacierwide_mb_mwea, elev_change_1d # class for Gaussian Process model for mass balance emulator @@ -726,6 +727,7 @@ def run(list_packed_vars): # ===== Load glacier data: area (km2), ice thickness (m), width (km) ===== try: + # Note this is where pre-processing of datasets (e.g., mass balance, debris) occurs if glacier_rgi_table['TermType'] not in [1, 5] or not pygem_prms['setup']['include_frontalablation']: gdir = oggm_compat.single_flowline_glacier_directory(glacier_str) gdir.is_tidewater = False @@ -1776,7 +1778,7 @@ def calc_mb_total_minelev(modelprms): ] # Melt # energy available for melt [degC day] - melt_energy_available = T_minelev * dates_table['daysinmonth'].values + melt_energy_available = T_minelev * dates_table['days_in_step'].values melt_energy_available[melt_energy_available < 0] = 0 # assume all snow melt because anything more would melt underlying ice in lowermost bin # SNOW MELT [m w.e.] diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py index aa252914..3a3bc0fd 100755 --- a/pygem/bin/run/run_simulation.py +++ b/pygem/bin/run/run_simulation.py @@ -17,7 +17,6 @@ # Built-in libraries import argparse import copy -import inspect import json import multiprocessing import os @@ -406,6 +405,10 @@ def run(list_packed_vars): print('sim years:', args.sim_startyear, args.sim_endyear) # ===== LOAD CLIMATE DATA ===== + # Future simulations are not yet set up with daily data + if pygem_prms['time']['timestep'] == 'daily': + assert args.sim_endyear <= 2025, 'Future daily data is not yet available' + # Climate class if sim_climate_name in ['ERA5', 'COAWST']: gcm = class_climate.GCM(name=sim_climate_name) @@ -482,6 +485,7 @@ def run(list_packed_vars): args.sim_startyear, args.ref_startyear, ) + # OPTION 2: Adjust temp and prec using Huss and Hock (2015) elif args.option_bias_adjustment == 2: # Temperature bias correction @@ -536,11 +540,13 @@ def run(list_packed_vars): gcm_tempstd = np.zeros((main_glac_rgi.shape[0], dates_table.shape[0])) ref_tempstd = np.zeros((main_glac_rgi.shape[0], dates_table_ref.shape[0])) elif pygem_prms['mb']['option_ablation'] == 2 and sim_climate_name in ['ERA5']: + assert pygem_prms['time']['timestep'] != 'daily', 'Option 2 for ablation should not be used with daily data' gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table, verbose=debug ) ref_tempstd = gcm_tempstd elif pygem_prms['mb']['option_ablation'] == 2 and args.ref_climate_name in ['ERA5']: + assert pygem_prms['time']['timestep'] != 'daily', 'Option 2 for ablation should not be used with daily data' # Compute temp std based on reference climate data ref_tempstd, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray( ref_gcm.tempstd_fn, @@ -556,29 +562,28 @@ def run(list_packed_vars): ref_tempstd = np.zeros((main_glac_rgi.shape[0], dates_table_ref.shape[0])) # Lapse rate - if sim_climate_name == 'ERA5': - gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( - gcm.lr_fn, - gcm.lr_vn, - main_glac_rgi, - dates_table, - upscale_var_timestep=True, - verbose=debug, - ) - ref_lr = gcm_lr + if pygem_prms['sim']['params']['use_constant_lapserate']: + gcm_lr = np.zeros((main_glac_rgi.shape[0], dates_table.shape[0])) + pygem_prms['sim']['params']['lapserate'] + ref_lr = np.zeros((main_glac_rgi.shape[0], dates_table_ref.shape[0])) + pygem_prms['sim']['params']['lapserate'] else: - # Compute lapse rates based on reference climate data - ref_lr, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray( - ref_gcm.lr_fn, ref_gcm.lr_vn, main_glac_rgi, dates_table_ref, verbose=debug - ) - # Monthly average from reference climate data - gcm_lr = gcmbiasadj.monthly_avg_array_rolled( - ref_lr, - dates_table_ref, - dates_table_full, - args.sim_startyear, - args.ref_startyear, - ) + if sim_climate_name in ['ERA-Interim', 'ERA5']: + gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( + gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table, upscale_var_timestep=True, verbose=debug + ) + ref_lr = gcm_lr + else: + # Compute lapse rates based on reference climate data + ref_lr, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray( + ref_gcm.lr_fn, ref_gcm.lr_vn, main_glac_rgi, dates_table_ref + ) + # Monthly average from reference climate data + gcm_lr = gcmbiasadj.monthly_avg_array_rolled( + ref_lr, + dates_table_ref, + dates_table_full, + args.sim_startyear, + args.ref_startyear, + ) # ===== RUN MASS BALANCE ===== # Number of simulations @@ -587,9 +592,9 @@ def run(list_packed_vars): else: nsims = 1 - # Number of years - nyears = dates_table.year.unique()[-1] - dates_table.year.unique()[0] + 1 - nyears_ref = dates_table_ref.year.unique()[-1] - dates_table.year.unique()[0] + 1 + # Number of years (for OGGM's run_until_and_store) + nyears = len(dates_table.year.unique()) + nyears_ref = len(dates_table_ref.year.unique()) for glac in range(main_glac_rgi.shape[0]): if glac == 0: @@ -605,8 +610,6 @@ def run(list_packed_vars): rgiid = main_glac_rgi.loc[main_glac_rgi.index.values[glac], 'RGIId'] try: - # for batman in [0]: - # ===== Load glacier data: area (km2), ice thickness (m), width (km) ===== if glacier_rgi_table['TermType'] not in [1, 5] or not pygem_prms['setup']['include_frontalablation']: gdir = single_flowline_glacier_directory(glacier_str, working_dir=args.oggm_working_dir) @@ -813,31 +816,32 @@ def run(list_packed_vars): # Time attributes and values if pygem_prms['climate']['sim_wateryear'] == 'hydro': - annual_columns = np.unique(dates_table['wateryear'].values)[0 : int(dates_table.shape[0] / 12)] + annual_columns = np.unique(dates_table['wateryear'].values) else: - annual_columns = np.unique(dates_table['year'].values)[0 : int(dates_table.shape[0] / 12)] + annual_columns = np.unique(dates_table['year'].values) + # append additional year to year_values to account for mass and area at end of period year_values = annual_columns year_values = np.concatenate((year_values, np.array([annual_columns[-1] + 1]))) - output_glac_temp_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_glac_prec_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_glac_acc_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_glac_refreeze_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_glac_melt_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_glac_frontalablation_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_glac_massbaltotal_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_glac_runoff_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_glac_snowline_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_glac_temp_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_glac_prec_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_glac_acc_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_glac_refreeze_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_glac_melt_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_glac_frontalablation_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_glac_massbaltotal_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_glac_runoff_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_glac_snowline_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan output_glac_area_annual = np.zeros((year_values.shape[0], nsims)) * np.nan output_glac_mass_annual = np.zeros((year_values.shape[0], nsims)) * np.nan output_glac_mass_bsl_annual = np.zeros((year_values.shape[0], nsims)) * np.nan output_glac_mass_change_ignored_annual = np.zeros((year_values.shape[0], nsims)) output_glac_ELA_annual = np.zeros((year_values.shape[0], nsims)) * np.nan - output_offglac_prec_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_offglac_refreeze_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_offglac_melt_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_offglac_snowpack_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan - output_offglac_runoff_monthly = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_offglac_prec_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_offglac_refreeze_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_offglac_melt_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_offglac_snowpack_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan + output_offglac_runoff_steps = np.zeros((dates_table.shape[0], nsims)) * np.nan output_glac_bin_icethickness_annual = None # Loop through model parameters @@ -875,7 +879,6 @@ def run(list_packed_vars): + str(np.round(modelprms['tbias'], 2)) ) - # %% # ----- ICE THICKNESS INVERSION using OGGM ----- if args.option_dynamics is not None: # Apply inversion_filter on mass balance with debris to avoid negative flux @@ -960,7 +963,9 @@ def run(list_packed_vars): option_areaconstant=False, ) - # Glacier dynamics model + ###################################### + ### OGGM dynamical evolution model ### + ###################################### if args.option_dynamics == 'OGGM': if debug: print('OGGM GLACIER DYNAMICS!') @@ -988,188 +993,64 @@ def run(list_packed_vars): if debug: fig, ax = plt.subplots(1) - graphics.plot_modeloutput_section(ev_model, ax=ax) - - try: - diag = ev_model.run_until_and_store(args.sim_endyear + 1) - ev_model.mb_model.glac_wide_volume_annual[-1] = diag.volume_m3[-1] - ev_model.mb_model.glac_wide_area_annual[-1] = diag.area_m2[-1] - - # Record frontal ablation for tidewater glaciers and update total mass balance - if gdir.is_tidewater: - # Glacier-wide frontal ablation (m3 w.e.) - # - note: diag.calving_m3 is cumulative calving - if debug: - print('\n\ndiag.calving_m3:', diag.calving_m3.values) - print( - 'calving_m3_since_y0:', - ev_model.calving_m3_since_y0, - ) - calving_m3_annual = ( - (diag.calving_m3.values[1:] - diag.calving_m3.values[0:-1]) - * pygem_prms['constants']['density_ice'] - / pygem_prms['constants']['density_water'] - ) - for n in np.arange(calving_m3_annual.shape[0]): - ev_model.mb_model.glac_wide_frontalablation[12 * n + 11] = calving_m3_annual[n] + graphics.plot_modeloutput_section( + ev_model, ax=ax, lnlabel=f'Glacier year {args.sim_startyear}' + ) - # Glacier-wide total mass balance (m3 w.e.) - ev_model.mb_model.glac_wide_massbaltotal = ( - ev_model.mb_model.glac_wide_massbaltotal - - ev_model.mb_model.glac_wide_frontalablation - ) + diag = ev_model.run_until_and_store(args.sim_endyear + 1) + ev_model.mb_model.glac_wide_volume_annual[-1] = diag.volume_m3[-1] + ev_model.mb_model.glac_wide_area_annual[-1] = diag.area_m2[-1] - if debug: - print( - 'avg calving_m3:', - calving_m3_annual.sum() / nyears, - ) - print( - 'avg frontal ablation [Gta]:', - np.round( - ev_model.mb_model.glac_wide_frontalablation.sum() / 1e9 / nyears, - 4, - ), - ) - print( - 'avg frontal ablation [Gta]:', - np.round( - ev_model.calving_m3_since_y0 - * pygem_prms['constants']['density_ice'] - / 1e12 - / nyears, - 4, - ), - ) - - except RuntimeError as e: - if 'Glacier exceeds domain boundaries' in repr(e): - count_exceed_boundary_errors += 1 - successful_run = False - - # LOG FAILURE - fail_domain_fp = ( - pygem_prms['root'] - + '/Output/simulations/fail-exceed_domain/' - + reg_str - + '/' - + sim_climate_name - + '/' - ) - if sim_climate_name not in [ - 'ERA5', - 'COAWST', - ]: - fail_domain_fp += sim_climate_scenario + '/' - if not os.path.exists(fail_domain_fp): - os.makedirs(fail_domain_fp, exist_ok=True) - txt_fn_fail = glacier_str + '-sim_failed.txt' - with open(fail_domain_fp + txt_fn_fail, 'w') as text_file: - text_file.write( - glacier_str - + ' failed to complete ' - + str(count_exceed_boundary_errors) - + ' simulations' - ) - elif gdir.is_tidewater: - if debug: - print('OGGM dynamics failed, using mass redistribution curves') - # Mass redistribution curves glacier dynamics model - ev_model = MassRedistributionCurveModel( - nfls, - mb_model=mbmod, - y0=args.sim_startyear, - glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier, - fs=fs, - is_tidewater=gdir.is_tidewater, - water_level=water_level, - ) - _, diag = ev_model.run_until_and_store(args.sim_endyear + 1) - ev_model.mb_model.glac_wide_volume_annual = diag.volume_m3.values - ev_model.mb_model.glac_wide_area_annual = diag.area_m2.values - - # Record frontal ablation for tidewater glaciers and update total mass balance - # Update glacier-wide frontal ablation (m3 w.e.) - ev_model.mb_model.glac_wide_frontalablation = ( - ev_model.mb_model.glac_bin_frontalablation.sum(0) - ) - # Update glacier-wide total mass balance (m3 w.e.) - ev_model.mb_model.glac_wide_massbaltotal = ( - ev_model.mb_model.glac_wide_massbaltotal - - ev_model.mb_model.glac_wide_frontalablation + # Record frontal ablation for tidewater glaciers and update total mass balance + if gdir.is_tidewater: + # Glacier-wide frontal ablation (m3 w.e.) + # - note: diag.calving_m3 is cumulative calving + if debug: + print('\n\ndiag.calving_m3:', diag.calving_m3.values) + print( + 'calving_m3_since_y0:', + ev_model.calving_m3_since_y0, ) + calving_m3_annual = ( + (diag.calving_m3.values[1:] - diag.calving_m3.values[0:-1]) + * pygem_prms['constants']['density_ice'] + / pygem_prms['constants']['density_water'] + ) + for n, year in enumerate(np.arange(args.sim_startyear, args.sim_endyear + 1)): + tstart, tstop = ev_model.mb_model.get_step_inds(year) + ev_model.mb_model.glac_wide_frontalablation[tstop] = calving_m3_annual[n] - if debug: - print( - 'avg frontal ablation [Gta]:', - np.round( - ev_model.mb_model.glac_wide_frontalablation.sum() / 1e9 / nyears, - 4, - ), - ) - print( - 'avg frontal ablation [Gta]:', - np.round( - ev_model.calving_m3_since_y0 - * pygem_prms['constants']['density_ice'] - / 1e12 - / nyears, - 4, - ), - ) + # Glacier-wide total mass balance (m3 w.e.) + ev_model.mb_model.glac_wide_massbaltotal = ( + ev_model.mb_model.glac_wide_massbaltotal - ev_model.mb_model.glac_wide_frontalablation + ) - except: - if gdir.is_tidewater: - if debug: - print('OGGM dynamics failed, using mass redistribution curves') - # Mass redistribution curves glacier dynamics model - ev_model = MassRedistributionCurveModel( - nfls, - mb_model=mbmod, - y0=args.sim_startyear, - glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier, - fs=fs, - is_tidewater=gdir.is_tidewater, - water_level=water_level, + if debug: + print( + 'avg calving_m3:', + calving_m3_annual.sum() / nyears, ) - _, diag = ev_model.run_until_and_store(args.sim_endyear + 1) - ev_model.mb_model.glac_wide_volume_annual = diag.volume_m3.values - ev_model.mb_model.glac_wide_area_annual = diag.area_m2.values - - # Record frontal ablation for tidewater glaciers and update total mass balance - # Update glacier-wide frontal ablation (m3 w.e.) - ev_model.mb_model.glac_wide_frontalablation = ( - ev_model.mb_model.glac_bin_frontalablation.sum(0) + print( + 'avg frontal ablation [Gta]:', + np.round( + ev_model.mb_model.glac_wide_frontalablation.sum() / 1e9 / nyears, + 4, + ), ) - # Update glacier-wide total mass balance (m3 w.e.) - ev_model.mb_model.glac_wide_massbaltotal = ( - ev_model.mb_model.glac_wide_massbaltotal - - ev_model.mb_model.glac_wide_frontalablation + print( + 'avg frontal ablation [Gta]:', + np.round( + ev_model.calving_m3_since_y0 + * pygem_prms['constants']['density_ice'] + / 1e12 + / nyears, + 4, + ), ) - if debug: - print( - 'avg frontal ablation [Gta]:', - np.round( - ev_model.mb_model.glac_wide_frontalablation.sum() / 1e9 / nyears, - 4, - ), - ) - print( - 'avg frontal ablation [Gta]:', - np.round( - ev_model.calving_m3_since_y0 - * pygem_prms['constants']['density_ice'] - / 1e12 - / nyears, - 4, - ), - ) - - else: - raise - - # Mass redistribution model + ###################################### + ##### mass redistribution model ##### + ###################################### elif args.option_dynamics == 'MassRedistributionCurves': if debug: print('MASS REDISTRIBUTION CURVES!') @@ -1180,82 +1061,54 @@ def run(list_packed_vars): glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier, fs=fs, is_tidewater=gdir.is_tidewater, - # water_level=gdir.get_diagnostics().get('calving_water_level', None) + # water_level=gdir.get_diagnostics().get('calving_water_level', None) water_level=water_level, ) if debug: fig, ax = plt.subplots(1) - graphics.plot_modeloutput_section(ev_model, ax=ax) - try: - _, diag = ev_model.run_until_and_store(args.sim_endyear + 1) - # print('shape of volume:', ev_model.mb_model.glac_wide_volume_annual.shape, diag.volume_m3.shape) - ev_model.mb_model.glac_wide_volume_annual = diag.volume_m3.values - ev_model.mb_model.glac_wide_area_annual = diag.area_m2.values - - # Record frontal ablation for tidewater glaciers and update total mass balance - if gdir.is_tidewater: - # Update glacier-wide frontal ablation (m3 w.e.) - ev_model.mb_model.glac_wide_frontalablation = ( - ev_model.mb_model.glac_bin_frontalablation.sum(0) - ) - # Update glacier-wide total mass balance (m3 w.e.) - ev_model.mb_model.glac_wide_massbaltotal = ( - ev_model.mb_model.glac_wide_massbaltotal - - ev_model.mb_model.glac_wide_frontalablation - ) + graphics.plot_modeloutput_section( + ev_model, ax=ax, lnlabel=f'Glacier year {args.sim_startyear}' + ) - if debug: - print( - 'avg frontal ablation [Gta]:', - np.round( - ev_model.mb_model.glac_wide_frontalablation.sum() / 1e9 / nyears, - 4, - ), - ) - print( - 'avg frontal ablation [Gta]:', - np.round( - ev_model.calving_m3_since_y0 - * pygem_prms['constants']['density_ice'] - / 1e12 - / nyears, - 4, - ), - ) + _, diag = ev_model.run_until_and_store(args.sim_endyear + 1) + # print('shape of volume:', ev_model.mb_model.glac_wide_volume_annual.shape, diag.volume_m3.shape) + ev_model.mb_model.glac_wide_volume_annual = diag.volume_m3.values + ev_model.mb_model.glac_wide_area_annual = diag.area_m2.values + + # Record frontal ablation for tidewater glaciers and update total mass balance + if gdir.is_tidewater: + # Update glacier-wide frontal ablation (m3 w.e.) + ev_model.mb_model.glac_wide_frontalablation = ( + ev_model.mb_model.glac_bin_frontalablation.sum(0) + ) + # Update glacier-wide total mass balance (m3 w.e.) + ev_model.mb_model.glac_wide_massbaltotal = ( + ev_model.mb_model.glac_wide_massbaltotal - ev_model.mb_model.glac_wide_frontalablation + ) - except RuntimeError as e: - if 'Glacier exceeds domain boundaries' in repr(e): - count_exceed_boundary_errors += 1 - successful_run = False - - # LOG FAILURE - fail_domain_fp = ( - pygem_prms['root'] - + '/Output/simulations/fail-exceed_domain/' - + reg_str - + '/' - + sim_climate_name - + '/' + if debug: + print( + 'avg frontal ablation [Gta]:', + np.round( + ev_model.mb_model.glac_wide_frontalablation.sum() / 1e9 / nyears, + 4, + ), + ) + print( + 'avg frontal ablation [Gta]:', + np.round( + ev_model.calving_m3_since_y0 + * pygem_prms['constants']['density_ice'] + / 1e12 + / nyears, + 4, + ), ) - if sim_climate_name not in [ - 'ERA5', - 'COAWST', - ]: - fail_domain_fp += sim_climate_scenario + '/' - if not os.path.exists(fail_domain_fp): - os.makedirs(fail_domain_fp, exist_ok=True) - txt_fn_fail = glacier_str + '-sim_failed.txt' - with open(fail_domain_fp + txt_fn_fail, 'w') as text_file: - text_file.write( - glacier_str - + ' failed to complete ' - + str(count_exceed_boundary_errors) - + ' simulations' - ) - else: - raise + ###################################### + ######### no dynamical model ######### + ###################################### elif args.option_dynamics is None: # Mass balance model ev_model = None @@ -1271,25 +1124,20 @@ def run(list_packed_vars): years = np.arange(args.sim_startyear, args.sim_endyear + 1) mb_all = [] for year in years: - mb_annual = mbmod.get_annual_mb( + # Calculate annual mass balance + mbmod.get_annual_mb( nfls[0].surface_h, fls=nfls, fl_id=0, year=year, - debug=True, + debug=debug, ) - mb_mwea = ( - mb_annual - * 365 - * 24 - * 3600 - * pygem_prms['constants']['density_ice'] - / pygem_prms['constants']['density_water'] + # Record glacierwide annual mass balance + t_start, t_stop = mbmod.get_step_inds(year) + mb_all.append( + mbmod.glac_wide_massbaltotal[t_start : t_stop + 1].sum() + / mbmod.glacier_area_initial.sum() ) - glac_wide_mb_mwea = ( - mb_mwea * mbmod.glacier_area_initial - ).sum() / mbmod.glacier_area_initial.sum() - mb_all.append(glac_wide_mb_mwea) mbmod.glac_wide_area_annual[-1] = mbmod.glac_wide_area_annual[0] mbmod.glac_wide_volume_annual[-1] = mbmod.glac_wide_volume_annual[0] diag['area_m2'] = mbmod.glac_wide_area_annual @@ -1307,15 +1155,13 @@ def run(list_packed_vars): np.round(np.median(mb_all), 3), ) - # mb_em_mwea = run_emulator_mb(modelprms) - # print(' emulator mb:', np.round(mb_em_mwea,3)) - # mb_em_sims.append(mb_em_mwea) - # Record output for successful runs if successful_run: if args.option_dynamics is not None: if debug: - graphics.plot_modeloutput_section(ev_model, ax=ax, srfls='--') + graphics.plot_modeloutput_section( + ev_model, ax=ax, srfls='--', lnlabel=f'Glacier year {args.sim_endyear + 1}' + ) plt.figure() diag.volume_m3.plot() plt.show() @@ -1359,15 +1205,15 @@ def run(list_packed_vars): ) # RECORD PARAMETERS TO DATASET - output_glac_temp_monthly[:, n_iter] = mbmod.glac_wide_temp - output_glac_prec_monthly[:, n_iter] = mbmod.glac_wide_prec - output_glac_acc_monthly[:, n_iter] = mbmod.glac_wide_acc - output_glac_refreeze_monthly[:, n_iter] = mbmod.glac_wide_refreeze - output_glac_melt_monthly[:, n_iter] = mbmod.glac_wide_melt - output_glac_frontalablation_monthly[:, n_iter] = mbmod.glac_wide_frontalablation - output_glac_massbaltotal_monthly[:, n_iter] = mbmod.glac_wide_massbaltotal - output_glac_runoff_monthly[:, n_iter] = mbmod.glac_wide_runoff - output_glac_snowline_monthly[:, n_iter] = mbmod.glac_wide_snowline + output_glac_temp_steps[:, n_iter] = mbmod.glac_wide_temp + output_glac_prec_steps[:, n_iter] = mbmod.glac_wide_prec + output_glac_acc_steps[:, n_iter] = mbmod.glac_wide_acc + output_glac_refreeze_steps[:, n_iter] = mbmod.glac_wide_refreeze + output_glac_melt_steps[:, n_iter] = mbmod.glac_wide_melt + output_glac_frontalablation_steps[:, n_iter] = mbmod.glac_wide_frontalablation + output_glac_massbaltotal_steps[:, n_iter] = mbmod.glac_wide_massbaltotal + output_glac_runoff_steps[:, n_iter] = mbmod.glac_wide_runoff + output_glac_snowline_steps[:, n_iter] = mbmod.glac_wide_snowline output_glac_area_annual[:, n_iter] = diag.area_m2.values output_glac_mass_annual[:, n_iter] = ( diag.volume_m3.values * pygem_prms['constants']['density_ice'] @@ -1379,12 +1225,12 @@ def run(list_packed_vars): mbmod.glac_wide_volume_change_ignored_annual * pygem_prms['constants']['density_ice'] ) output_glac_ELA_annual[:, n_iter] = mbmod.glac_wide_ELA_annual - output_offglac_prec_monthly[:, n_iter] = mbmod.offglac_wide_prec + output_offglac_prec_steps[:, n_iter] = mbmod.offglac_wide_prec - output_offglac_refreeze_monthly[:, n_iter] = mbmod.offglac_wide_refreeze - output_offglac_melt_monthly[:, n_iter] = mbmod.offglac_wide_melt - output_offglac_snowpack_monthly[:, n_iter] = mbmod.offglac_wide_snowpack - output_offglac_runoff_monthly[:, n_iter] = mbmod.offglac_wide_runoff + output_offglac_refreeze_steps[:, n_iter] = mbmod.offglac_wide_refreeze + output_offglac_melt_steps[:, n_iter] = mbmod.offglac_wide_melt + output_offglac_snowpack_steps[:, n_iter] = mbmod.offglac_wide_snowpack + output_offglac_runoff_steps[:, n_iter] = mbmod.offglac_wide_runoff if output_glac_bin_icethickness_annual is None: output_glac_bin_area_annual_sim = mbmod.glac_bin_area_annual[:, :, np.newaxis] @@ -1425,23 +1271,21 @@ def run(list_packed_vars): output_glac_bin_massbalclim_annual = output_glac_bin_massbalclim_annual_sim[ :, :, np.newaxis ] - output_glac_bin_massbalclim_monthly_sim = np.zeros(mbmod.glac_bin_massbalclim.shape) - output_glac_bin_massbalclim_monthly_sim = mbmod.glac_bin_massbalclim - output_glac_bin_massbalclim_monthly = output_glac_bin_massbalclim_monthly_sim[ - :, :, np.newaxis - ] - # accum - output_glac_bin_acc_monthly_sim = np.zeros(mbmod.bin_acc.shape) - output_glac_bin_acc_monthly_sim = mbmod.bin_acc - output_glac_bin_acc_monthly = output_glac_bin_acc_monthly_sim[:, :, np.newaxis] + output_glac_bin_massbalclim_steps_sim = np.zeros(mbmod.glac_bin_massbalclim.shape) + output_glac_bin_massbalclim_steps_sim = mbmod.glac_bin_massbalclim + output_glac_bin_massbalclim_steps = output_glac_bin_massbalclim_steps_sim[:, :, np.newaxis] + # accumulation + output_glac_bin_acc_steps_sim = np.zeros(mbmod.bin_acc.shape) + output_glac_bin_acc_steps_sim = mbmod.bin_acc + output_glac_bin_acc_steps = output_glac_bin_acc_steps_sim[:, :, np.newaxis] # refreeze - output_glac_bin_refreeze_monthly_sim = np.zeros(mbmod.glac_bin_refreeze.shape) - output_glac_bin_refreeze_monthly_sim = mbmod.glac_bin_refreeze - output_glac_bin_refreeze_monthly = output_glac_bin_refreeze_monthly_sim[:, :, np.newaxis] + output_glac_bin_refreeze_steps_sim = np.zeros(mbmod.glac_bin_refreeze.shape) + output_glac_bin_refreeze_steps_sim = mbmod.glac_bin_refreeze + output_glac_bin_refreeze_steps = output_glac_bin_refreeze_steps_sim[:, :, np.newaxis] # melt - output_glac_bin_melt_monthly_sim = np.zeros(mbmod.glac_bin_melt.shape) - output_glac_bin_melt_monthly_sim = mbmod.glac_bin_melt - output_glac_bin_melt_monthly = output_glac_bin_melt_monthly_sim[:, :, np.newaxis] + output_glac_bin_melt_steps_sim = np.zeros(mbmod.glac_bin_melt.shape) + output_glac_bin_melt_steps_sim = mbmod.glac_bin_melt + output_glac_bin_melt_steps = output_glac_bin_melt_steps_sim[:, :, np.newaxis] else: # Update the latest thickness and volume @@ -1496,35 +1340,35 @@ def run(list_packed_vars): output_glac_bin_massbalclim_annual_sim[:, :, np.newaxis], axis=2, ) - output_glac_bin_massbalclim_monthly_sim = np.zeros(mbmod.glac_bin_massbalclim.shape) - output_glac_bin_massbalclim_monthly_sim = mbmod.glac_bin_massbalclim - output_glac_bin_massbalclim_monthly = np.append( - output_glac_bin_massbalclim_monthly, - output_glac_bin_massbalclim_monthly_sim[:, :, np.newaxis], + output_glac_bin_massbalclim_steps_sim = np.zeros(mbmod.glac_bin_massbalclim.shape) + output_glac_bin_massbalclim_steps_sim = mbmod.glac_bin_massbalclim + output_glac_bin_massbalclim_steps = np.append( + output_glac_bin_massbalclim_steps, + output_glac_bin_massbalclim_steps_sim[:, :, np.newaxis], axis=2, ) - # accum - output_glac_bin_acc_monthly_sim = np.zeros(mbmod.bin_acc.shape) - output_glac_bin_acc_monthly_sim = mbmod.bin_acc - output_glac_bin_acc_monthly = np.append( - output_glac_bin_acc_monthly, - output_glac_bin_acc_monthly_sim[:, :, np.newaxis], + # accumulation + output_glac_bin_acc_steps_sim = np.zeros(mbmod.bin_acc.shape) + output_glac_bin_acc_steps_sim = mbmod.bin_acc + output_glac_bin_acc_steps = np.append( + output_glac_bin_acc_steps, + output_glac_bin_acc_steps_sim[:, :, np.newaxis], axis=2, ) # melt - output_glac_bin_melt_monthly_sim = np.zeros(mbmod.glac_bin_melt.shape) - output_glac_bin_melt_monthly_sim = mbmod.glac_bin_melt - output_glac_bin_melt_monthly = np.append( - output_glac_bin_melt_monthly, - output_glac_bin_melt_monthly_sim[:, :, np.newaxis], + output_glac_bin_melt_steps_sim = np.zeros(mbmod.glac_bin_melt.shape) + output_glac_bin_melt_steps_sim = mbmod.glac_bin_melt + output_glac_bin_melt_steps = np.append( + output_glac_bin_melt_steps, + output_glac_bin_melt_steps_sim[:, :, np.newaxis], axis=2, ) # refreeze - output_glac_bin_refreeze_monthly_sim = np.zeros(mbmod.glac_bin_refreeze.shape) - output_glac_bin_refreeze_monthly_sim = mbmod.glac_bin_refreeze - output_glac_bin_refreeze_monthly = np.append( - output_glac_bin_refreeze_monthly, - output_glac_bin_refreeze_monthly_sim[:, :, np.newaxis], + output_glac_bin_refreeze_steps_sim = np.zeros(mbmod.glac_bin_refreeze.shape) + output_glac_bin_refreeze_steps_sim = mbmod.glac_bin_refreeze + output_glac_bin_refreeze_steps = np.append( + output_glac_bin_refreeze_steps, + output_glac_bin_refreeze_steps_sim[:, :, np.newaxis], axis=2, ) @@ -1537,6 +1381,7 @@ def run(list_packed_vars): output_stats = output.glacierwide_stats( glacier_rgi_table=glacier_rgi_table, dates_table=dates_table, + timestep=pygem_prms['time']['timestep'], nsims=1, sim_climate_name=sim_climate_name, sim_climate_scenario=sim_climate_scenario, @@ -1548,8 +1393,12 @@ def run(list_packed_vars): sim_endyear=args.sim_endyear, option_calibration=args.option_calibration, option_bias_adjustment=args.option_bias_adjustment, + option_dynamics=args.option_dynamics, extra_vars=args.export_extra_vars, ) + base_fn = ( + output_stats.get_fn() + ) # should contain 'SETS' which is later used to replace with the specific iteration for n_iter in range(nsims): # pass model params for iteration and update output dataset model params output_stats.set_modelprms({key: modelprms_all[key][n_iter] for key in modelprms_all}) @@ -1557,69 +1406,54 @@ def run(list_packed_vars): output_stats.create_xr_ds() output_ds_all_stats = output_stats.get_xr_ds() # fill values - output_ds_all_stats['glac_runoff_monthly'].values[0, :] = output_glac_runoff_monthly[ - :, n_iter - ] + output_ds_all_stats['glac_runoff'].values[0, :] = output_glac_runoff_steps[:, n_iter] output_ds_all_stats['glac_area_annual'].values[0, :] = output_glac_area_annual[:, n_iter] output_ds_all_stats['glac_mass_annual'].values[0, :] = output_glac_mass_annual[:, n_iter] output_ds_all_stats['glac_mass_bsl_annual'].values[0, :] = output_glac_mass_bsl_annual[ :, n_iter ] output_ds_all_stats['glac_ELA_annual'].values[0, :] = output_glac_ELA_annual[:, n_iter] - output_ds_all_stats['offglac_runoff_monthly'].values[0, :] = output_offglac_runoff_monthly[ - :, n_iter - ] + output_ds_all_stats['offglac_runoff'].values[0, :] = output_offglac_runoff_steps[:, n_iter] if args.export_extra_vars: - output_ds_all_stats['glac_temp_monthly'].values[0, :] = ( - output_glac_temp_monthly[:, n_iter] + 273.15 + output_ds_all_stats['glac_temp'].values[0, :] = ( + output_glac_temp_steps[:, n_iter] + 273.15 ) - output_ds_all_stats['glac_prec_monthly'].values[0, :] = output_glac_prec_monthly[ + output_ds_all_stats['glac_prec'].values[0, :] = output_glac_prec_steps[:, n_iter] + output_ds_all_stats['glac_acc'].values[0, :] = output_glac_acc_steps[:, n_iter] + output_ds_all_stats['glac_refreeze'].values[0, :] = output_glac_refreeze_steps[ :, n_iter ] - output_ds_all_stats['glac_acc_monthly'].values[0, :] = output_glac_acc_monthly[ + output_ds_all_stats['glac_melt'].values[0, :] = output_glac_melt_steps[:, n_iter] + output_ds_all_stats['glac_frontalablation'].values[0, :] = ( + output_glac_frontalablation_steps[:, n_iter] + ) + output_ds_all_stats['glac_massbaltotal'].values[0, :] = output_glac_massbaltotal_steps[ :, n_iter ] - output_ds_all_stats['glac_refreeze_monthly'].values[0, :] = ( - output_glac_refreeze_monthly[:, n_iter] - ) - output_ds_all_stats['glac_melt_monthly'].values[0, :] = output_glac_melt_monthly[ + output_ds_all_stats['glac_snowline'].values[0, :] = output_glac_snowline_steps[ :, n_iter ] - output_ds_all_stats['glac_frontalablation_monthly'].values[0, :] = ( - output_glac_frontalablation_monthly[:, n_iter] - ) - output_ds_all_stats['glac_massbaltotal_monthly'].values[0, :] = ( - output_glac_massbaltotal_monthly[:, n_iter] - ) - output_ds_all_stats['glac_snowline_monthly'].values[0, :] = ( - output_glac_snowline_monthly[:, n_iter] - ) output_ds_all_stats['glac_mass_change_ignored_annual'].values[0, :] = ( output_glac_mass_change_ignored_annual[:, n_iter] ) - output_ds_all_stats['offglac_prec_monthly'].values[0, :] = output_offglac_prec_monthly[ + output_ds_all_stats['offglac_prec'].values[0, :] = output_offglac_prec_steps[:, n_iter] + output_ds_all_stats['offglac_melt'].values[0, :] = output_offglac_melt_steps[:, n_iter] + output_ds_all_stats['offglac_refreeze'].values[0, :] = output_offglac_refreeze_steps[ :, n_iter ] - output_ds_all_stats['offglac_melt_monthly'].values[0, :] = output_offglac_melt_monthly[ + output_ds_all_stats['offglac_snowpack'].values[0, :] = output_offglac_snowpack_steps[ :, n_iter ] - output_ds_all_stats['offglac_refreeze_monthly'].values[0, :] = ( - output_offglac_refreeze_monthly[:, n_iter] - ) - output_ds_all_stats['offglac_snowpack_monthly'].values[0, :] = ( - output_offglac_snowpack_monthly[:, n_iter] - ) # export glacierwide stats for iteration - output_stats.set_fn( - output_stats.get_fn().replace('SETS', f'set{n_iter}') + args.outputfn_sfix + 'all.nc' - ) + output_stats.set_fn(base_fn.replace('SETS', f'set{n_iter}') + args.outputfn_sfix + 'all.nc') output_stats.save_xr_ds() # instantiate dataset for merged simulations output_stats = output.glacierwide_stats( glacier_rgi_table=glacier_rgi_table, dates_table=dates_table, + timestep=pygem_prms['time']['timestep'], nsims=nsims, sim_climate_name=sim_climate_name, sim_climate_scenario=sim_climate_scenario, @@ -1631,6 +1465,7 @@ def run(list_packed_vars): sim_endyear=args.sim_endyear, option_calibration=args.option_calibration, option_bias_adjustment=args.option_bias_adjustment, + option_dynamics=args.option_dynamics, extra_vars=args.export_extra_vars, ) # create and return xarray dataset @@ -1638,128 +1473,95 @@ def run(list_packed_vars): output_ds_all_stats = output_stats.get_xr_ds() # get stats from all simulations which will be stored - output_glac_runoff_monthly_stats = calc_stats_array(output_glac_runoff_monthly) + output_glac_runoff_steps_stats = calc_stats_array(output_glac_runoff_steps) output_glac_area_annual_stats = calc_stats_array(output_glac_area_annual) output_glac_mass_annual_stats = calc_stats_array(output_glac_mass_annual) output_glac_mass_bsl_annual_stats = calc_stats_array(output_glac_mass_bsl_annual) output_glac_ELA_annual_stats = calc_stats_array(output_glac_ELA_annual) - output_offglac_runoff_monthly_stats = calc_stats_array(output_offglac_runoff_monthly) + output_offglac_runoff_steps_stats = calc_stats_array(output_offglac_runoff_steps) if args.export_extra_vars: - output_glac_temp_monthly_stats = calc_stats_array(output_glac_temp_monthly) - output_glac_prec_monthly_stats = calc_stats_array(output_glac_prec_monthly) - output_glac_acc_monthly_stats = calc_stats_array(output_glac_acc_monthly) - output_glac_refreeze_monthly_stats = calc_stats_array(output_glac_refreeze_monthly) - output_glac_melt_monthly_stats = calc_stats_array(output_glac_melt_monthly) - output_glac_frontalablation_monthly_stats = calc_stats_array( - output_glac_frontalablation_monthly - ) - output_glac_massbaltotal_monthly_stats = calc_stats_array(output_glac_massbaltotal_monthly) - output_glac_snowline_monthly_stats = calc_stats_array(output_glac_snowline_monthly) + output_glac_temp_steps_stats = calc_stats_array(output_glac_temp_steps) + output_glac_prec_steps_stats = calc_stats_array(output_glac_prec_steps) + output_glac_acc_steps_stats = calc_stats_array(output_glac_acc_steps) + output_glac_refreeze_steps_stats = calc_stats_array(output_glac_refreeze_steps) + output_glac_melt_steps_stats = calc_stats_array(output_glac_melt_steps) + output_glac_frontalablation_steps_stats = calc_stats_array(output_glac_frontalablation_steps) + output_glac_massbaltotal_steps_stats = calc_stats_array(output_glac_massbaltotal_steps) + output_glac_snowline_steps_stats = calc_stats_array(output_glac_snowline_steps) output_glac_mass_change_ignored_annual_stats = calc_stats_array( output_glac_mass_change_ignored_annual ) - output_offglac_prec_monthly_stats = calc_stats_array(output_offglac_prec_monthly) - output_offglac_melt_monthly_stats = calc_stats_array(output_offglac_melt_monthly) - output_offglac_refreeze_monthly_stats = calc_stats_array(output_offglac_refreeze_monthly) - output_offglac_snowpack_monthly_stats = calc_stats_array(output_offglac_snowpack_monthly) + output_offglac_prec_steps_stats = calc_stats_array(output_offglac_prec_steps) + output_offglac_melt_steps_stats = calc_stats_array(output_offglac_melt_steps) + output_offglac_refreeze_steps_stats = calc_stats_array(output_offglac_refreeze_steps) + output_offglac_snowpack_steps_stats = calc_stats_array(output_offglac_snowpack_steps) # output mean/median from all simulations - output_ds_all_stats['glac_runoff_monthly'].values[0, :] = output_glac_runoff_monthly_stats[:, 0] + output_ds_all_stats['glac_runoff'].values[0, :] = output_glac_runoff_steps_stats[:, 0] output_ds_all_stats['glac_area_annual'].values[0, :] = output_glac_area_annual_stats[:, 0] output_ds_all_stats['glac_mass_annual'].values[0, :] = output_glac_mass_annual_stats[:, 0] output_ds_all_stats['glac_mass_bsl_annual'].values[0, :] = output_glac_mass_bsl_annual_stats[:, 0] output_ds_all_stats['glac_ELA_annual'].values[0, :] = output_glac_ELA_annual_stats[:, 0] - output_ds_all_stats['offglac_runoff_monthly'].values[0, :] = output_offglac_runoff_monthly_stats[ - :, 0 - ] + output_ds_all_stats['offglac_runoff'].values[0, :] = output_offglac_runoff_steps_stats[:, 0] if args.export_extra_vars: - output_ds_all_stats['glac_temp_monthly'].values[0, :] = ( - output_glac_temp_monthly_stats[:, 0] + 273.15 - ) - output_ds_all_stats['glac_prec_monthly'].values[0, :] = output_glac_prec_monthly_stats[:, 0] - output_ds_all_stats['glac_acc_monthly'].values[0, :] = output_glac_acc_monthly_stats[:, 0] - output_ds_all_stats['glac_refreeze_monthly'].values[0, :] = output_glac_refreeze_monthly_stats[ - :, 0 - ] - output_ds_all_stats['glac_melt_monthly'].values[0, :] = output_glac_melt_monthly_stats[:, 0] - output_ds_all_stats['glac_frontalablation_monthly'].values[0, :] = ( - output_glac_frontalablation_monthly_stats[:, 0] + output_ds_all_stats['glac_temp'].values[0, :] = output_glac_temp_steps_stats[:, 0] + 273.15 + output_ds_all_stats['glac_prec'].values[0, :] = output_glac_prec_steps_stats[:, 0] + output_ds_all_stats['glac_acc'].values[0, :] = output_glac_acc_steps_stats[:, 0] + output_ds_all_stats['glac_refreeze'].values[0, :] = output_glac_refreeze_steps_stats[:, 0] + output_ds_all_stats['glac_melt'].values[0, :] = output_glac_melt_steps_stats[:, 0] + output_ds_all_stats['glac_frontalablation'].values[0, :] = ( + output_glac_frontalablation_steps_stats[:, 0] ) - output_ds_all_stats['glac_massbaltotal_monthly'].values[0, :] = ( - output_glac_massbaltotal_monthly_stats[:, 0] - ) - output_ds_all_stats['glac_snowline_monthly'].values[0, :] = output_glac_snowline_monthly_stats[ + output_ds_all_stats['glac_massbaltotal'].values[0, :] = output_glac_massbaltotal_steps_stats[ :, 0 ] + output_ds_all_stats['glac_snowline'].values[0, :] = output_glac_snowline_steps_stats[:, 0] output_ds_all_stats['glac_mass_change_ignored_annual'].values[0, :] = ( output_glac_mass_change_ignored_annual_stats[:, 0] ) - output_ds_all_stats['offglac_prec_monthly'].values[0, :] = output_offglac_prec_monthly_stats[ - :, 0 - ] - output_ds_all_stats['offglac_melt_monthly'].values[0, :] = output_offglac_melt_monthly_stats[ - :, 0 - ] - output_ds_all_stats['offglac_refreeze_monthly'].values[0, :] = ( - output_offglac_refreeze_monthly_stats[:, 0] - ) - output_ds_all_stats['offglac_snowpack_monthly'].values[0, :] = ( - output_offglac_snowpack_monthly_stats[:, 0] - ) + output_ds_all_stats['offglac_prec'].values[0, :] = output_offglac_prec_steps_stats[:, 0] + output_ds_all_stats['offglac_melt'].values[0, :] = output_offglac_melt_steps_stats[:, 0] + output_ds_all_stats['offglac_refreeze'].values[0, :] = output_offglac_refreeze_steps_stats[:, 0] + output_ds_all_stats['offglac_snowpack'].values[0, :] = output_offglac_snowpack_steps_stats[:, 0] # output median absolute deviation if nsims > 1: - output_ds_all_stats['glac_runoff_monthly_mad'].values[0, :] = output_glac_runoff_monthly_stats[ - :, 1 - ] + output_ds_all_stats['glac_runoff_mad'].values[0, :] = output_glac_runoff_steps_stats[:, 1] output_ds_all_stats['glac_area_annual_mad'].values[0, :] = output_glac_area_annual_stats[:, 1] output_ds_all_stats['glac_mass_annual_mad'].values[0, :] = output_glac_mass_annual_stats[:, 1] output_ds_all_stats['glac_mass_bsl_annual_mad'].values[0, :] = ( output_glac_mass_bsl_annual_stats[:, 1] ) output_ds_all_stats['glac_ELA_annual_mad'].values[0, :] = output_glac_ELA_annual_stats[:, 1] - output_ds_all_stats['offglac_runoff_monthly_mad'].values[0, :] = ( - output_offglac_runoff_monthly_stats[:, 1] - ) + output_ds_all_stats['offglac_runoff_mad'].values[0, :] = output_offglac_runoff_steps_stats[:, 1] + output_ds_all_stats['offglac_runoff_mad'].values[0, :] = output_offglac_runoff_steps_stats[:, 1] if args.export_extra_vars: - output_ds_all_stats['glac_temp_monthly_mad'].values[0, :] = output_glac_temp_monthly_stats[ + output_ds_all_stats['glac_temp_mad'].values[0, :] = output_glac_temp_steps_stats[:, 1] + output_ds_all_stats['glac_prec_mad'].values[0, :] = output_glac_prec_steps_stats[:, 1] + output_ds_all_stats['glac_acc_mad'].values[0, :] = output_glac_acc_steps_stats[:, 1] + output_ds_all_stats['glac_refreeze_mad'].values[0, :] = output_glac_refreeze_steps_stats[ :, 1 ] - output_ds_all_stats['glac_prec_monthly_mad'].values[0, :] = output_glac_prec_monthly_stats[ - :, 1 - ] - output_ds_all_stats['glac_acc_monthly_mad'].values[0, :] = output_glac_acc_monthly_stats[ - :, 1 - ] - output_ds_all_stats['glac_refreeze_monthly_mad'].values[0, :] = ( - output_glac_refreeze_monthly_stats[:, 1] + output_ds_all_stats['glac_melt_mad'].values[0, :] = output_glac_melt_steps_stats[:, 1] + output_ds_all_stats['glac_frontalablation_mad'].values[0, :] = ( + output_glac_frontalablation_steps_stats[:, 1] ) - output_ds_all_stats['glac_melt_monthly_mad'].values[0, :] = output_glac_melt_monthly_stats[ + output_ds_all_stats['glac_massbaltotal_mad'].values[0, :] = ( + output_glac_massbaltotal_steps_stats[:, 1] + ) + output_ds_all_stats['glac_snowline_mad'].values[0, :] = output_glac_snowline_steps_stats[ :, 1 ] - output_ds_all_stats['glac_frontalablation_monthly_mad'].values[0, :] = ( - output_glac_frontalablation_monthly_stats[:, 1] - ) - output_ds_all_stats['glac_massbaltotal_monthly_mad'].values[0, :] = ( - output_glac_massbaltotal_monthly_stats[:, 1] - ) - output_ds_all_stats['glac_snowline_monthly_mad'].values[0, :] = ( - output_glac_snowline_monthly_stats[:, 1] - ) output_ds_all_stats['glac_mass_change_ignored_annual_mad'].values[0, :] = ( output_glac_mass_change_ignored_annual_stats[:, 1] ) - output_ds_all_stats['offglac_prec_monthly_mad'].values[0, :] = ( - output_offglac_prec_monthly_stats[:, 1] - ) - output_ds_all_stats['offglac_melt_monthly_mad'].values[0, :] = ( - output_offglac_melt_monthly_stats[:, 1] + output_ds_all_stats['offglac_prec_mad'].values[0, :] = output_offglac_prec_steps_stats[:, 1] + output_ds_all_stats['offglac_melt_mad'].values[0, :] = output_offglac_melt_steps_stats[:, 1] + output_ds_all_stats['offglac_refreeze_mad'].values[0, :] = ( + output_offglac_refreeze_steps_stats[:, 1] ) - output_ds_all_stats['offglac_refreeze_monthly_mad'].values[0, :] = ( - output_offglac_refreeze_monthly_stats[:, 1] - ) - output_ds_all_stats['offglac_snowpack_monthly_mad'].values[0, :] = ( - output_offglac_snowpack_monthly_stats[:, 1] + output_ds_all_stats['offglac_snowpack_mad'].values[0, :] = ( + output_offglac_snowpack_steps_stats[:, 1] ) # export merged netcdf glacierwide stats @@ -1781,6 +1583,7 @@ def run(list_packed_vars): output_binned = output.binned_stats( glacier_rgi_table=glacier_rgi_table, dates_table=dates_table, + timestep=pygem_prms['time']['timestep'], nsims=1, nbins=surface_h_initial.shape[0], binned_components=args.export_binned_components, @@ -1794,7 +1597,11 @@ def run(list_packed_vars): sim_endyear=args.sim_endyear, option_calibration=args.option_calibration, option_bias_adjustment=args.option_bias_adjustment, + option_dynamics=args.option_dynamics, ) + base_fn = ( + output_binned.get_fn() + ) # should contain 'SETS' which is later used to replace with the specific iteration for n_iter in range(nsims): # pass model params for iteration and update output dataset model params output_binned.set_modelprms({key: modelprms_all[key][n_iter] for key in modelprms_all}) @@ -1816,25 +1623,23 @@ def run(list_packed_vars): output_ds_binned_stats['bin_massbalclim_annual'].values[0, :, :] = ( output_glac_bin_massbalclim_annual[:, :, n_iter] ) - output_ds_binned_stats['bin_massbalclim_monthly'].values[0, :, :] = ( - output_glac_bin_massbalclim_monthly[:, :, n_iter] + output_ds_binned_stats['bin_massbalclim'].values[0, :, :] = ( + output_glac_bin_massbalclim_steps[:, :, n_iter] ) if args.export_binned_components: - output_ds_binned_stats['bin_accumulation_monthly'].values[0, :, :] = ( - output_glac_bin_acc_monthly[:, :, n_iter] - ) - output_ds_binned_stats['bin_melt_monthly'].values[0, :, :] = ( - output_glac_bin_melt_monthly[:, :, n_iter] + output_ds_binned_stats['bin_accumulation'].values[0, :, :] = ( + output_glac_bin_acc_steps[:, :, n_iter] ) - output_ds_binned_stats['bin_refreeze_monthly'].values[0, :, :] = ( - output_glac_bin_refreeze_monthly[:, :, n_iter] + output_ds_binned_stats['bin_melt'].values[0, :, :] = output_glac_bin_melt_steps[ + :, :, n_iter + ] + output_ds_binned_stats['bin_refreeze'].values[0, :, :] = ( + output_glac_bin_refreeze_steps[:, :, n_iter] ) # export binned stats for iteration output_binned.set_fn( - output_binned.get_fn().replace('SETS', f'set{n_iter}') - + args.outputfn_sfix - + 'binned.nc' + base_fn.replace('SETS', f'set{n_iter}') + args.outputfn_sfix + 'binned.nc' ) output_binned.save_xr_ds() @@ -1842,6 +1647,7 @@ def run(list_packed_vars): output_binned = output.binned_stats( glacier_rgi_table=glacier_rgi_table, dates_table=dates_table, + timestep=pygem_prms['time']['timestep'], nsims=nsims, nbins=surface_h_initial.shape[0], binned_components=args.export_binned_components, @@ -1855,6 +1661,7 @@ def run(list_packed_vars): sim_endyear=args.sim_endyear, option_calibration=args.option_calibration, option_bias_adjustment=args.option_bias_adjustment, + option_dynamics=args.option_dynamics, ) # create and return xarray dataset output_binned.create_xr_ds() @@ -1875,18 +1682,18 @@ def run(list_packed_vars): output_ds_binned_stats['bin_massbalclim_annual'].values = np.median( output_glac_bin_massbalclim_annual, axis=2 )[np.newaxis, :, :] - output_ds_binned_stats['bin_massbalclim_monthly'].values = np.median( - output_glac_bin_massbalclim_monthly, axis=2 + output_ds_binned_stats['bin_massbalclim'].values = np.median( + output_glac_bin_massbalclim_steps, axis=2 )[np.newaxis, :, :] if args.export_binned_components: - output_ds_binned_stats['bin_accumulation_monthly'].values = np.median( - output_glac_bin_acc_monthly, axis=2 + output_ds_binned_stats['bin_accumulation'].values = np.median( + output_glac_bin_acc_steps, axis=2 )[np.newaxis, :, :] - output_ds_binned_stats['bin_melt_monthly'].values = np.median( - output_glac_bin_melt_monthly, axis=2 - )[np.newaxis, :, :] - output_ds_binned_stats['bin_refreeze_monthly'].values = np.median( - output_glac_bin_refreeze_monthly, axis=2 + output_ds_binned_stats['bin_melt'].values = np.median(output_glac_bin_melt_steps, axis=2)[ + np.newaxis, :, : + ] + output_ds_binned_stats['bin_refreeze'].values = np.median( + output_glac_bin_refreeze_steps, axis=2 )[np.newaxis, :, :] if nsims > 1: output_ds_binned_stats['bin_mass_annual_mad'].values = median_abs_deviation( @@ -1916,11 +1723,6 @@ def run(list_packed_vars): with open(fail_fp + txt_fn_fail, 'w') as text_file: text_file.write(glacier_str + f' failed to complete simulation: {err}') - # Global variables for Spyder development - if args.ncores == 1: - global main_vars - main_vars = inspect.currentframe().f_locals - # %% PARALLEL PROCESSING def main(): diff --git a/pygem/bin/run/run_spinup.py b/pygem/bin/run/run_spinup.py index 9d3d03f3..3b0e9fa5 100644 --- a/pygem/bin/run/run_spinup.py +++ b/pygem/bin/run/run_spinup.py @@ -34,8 +34,11 @@ from pygem.utils._funcs import interp1d_fill_gaps -# get model monthly deltah -def get_elev_change_1d_hat(gdir): +def calculate_elev_change_1d(gdir): + """ + calculate the 1d elevation change from dynamical spinup. + assumes constant flux divergence througohut each model year. + """ # load flowline_diagnostics from spinup f = gdir.get_filepath('fl_diagnostics', filesuffix='_dynamic_spinup_pygem_mb') with xr.open_dataset(f, group='fl_0') as ds_spn: @@ -50,55 +53,65 @@ def get_elev_change_1d_hat(gdir): 0, ) - thickness_m = ds_spn.thickness_m.values.T # glacier thickness [m ice], (nbins, nyears) - - # set any < 0 thickness to nan - thickness_m[thickness_m <= 0] = np.nan - - # climatic mass balance - dotb_monthly = np.repeat(ds_spn.climatic_mb_myr.values.T[:, 1:] / 12, 12, axis=-1) - - # convert to m ice - dotb_monthly = dotb_monthly * (pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice']) - ### to get monthly thickness and mass we require monthly flux divergence ### - # we'll assume the flux divergence is constant througohut the year (is this a good assumption?) - # ie. take annual values and divide by 12 - use numpy repeat to repeat values across 12 months - flux_div_monthly_mmo = np.repeat(-ds_spn.flux_divergence_myr.values.T[:, 1:] / 12, 12, axis=-1) - # get monthly binned change in thickness - delta_h_monthly = dotb_monthly - flux_div_monthly_mmo # [m ice per month] - - # get binned monthly thickness = running thickness change + initial thickness - running_delta_h_monthly = np.cumsum(delta_h_monthly, axis=-1) - h_monthly = running_delta_h_monthly + thickness_m[:, 0][:, np.newaxis] - + ### get monthly elevation change ### + # note, transpose and [:, 1:] indexing to access desired OGGM dataset outputs + # OGGM stores data following (nsteps, nbins) - while PyGEM follows (nbins, nsteps), hence .T operator + # OGGM also stores the 0th model year with np.nan values, hence [:,1:] indexing. Time index 1 corresponds to the output from model year 0-1 + # --- Step 1: convert mass balance from m w.e. to m ice --- + rho_w = pygem_prms['constants']['density_water'] + rho_i = pygem_prms['constants']['density_ice'] + bin_massbalclim_ice = ds_spn.climatic_mb_myr.values.T[:, 1:] * ( + rho_w / rho_i + ) # binned climatic mass balance (nbins, nsteps) + + # --- Step 2: expand annual climatic mass balance and flux divergence to monthly steps --- + # assume the climatic mass balance and flux divergence are constant througohut the year + # ie. take annual values and divide spread uniformly throughout model year + bin_massbalclim_ice_monthly = np.repeat(bin_massbalclim_ice / 12, 12, axis=-1) # [m ice] + bin_flux_divergence_monthly = np.repeat( + -ds_spn.flux_divergence_myr.values.T[:, 1:] / 12, 12, axis=-1 + ) # [m ice] note, oggm flux_divergence_myr is opposite sign of convention, hence negative + + # --- Step 3: compute monthly thickness change --- + bin_delta_thick_monthly = bin_massbalclim_ice_monthly - bin_flux_divergence_monthly # [m ice] + + # --- Step 4: calculate monthly thickness --- + # calculate binned monthly thickness = running thickness change + initial thickness + bin_thick_initial = ds_spn.thickness_m.isel(time=0).values # initial glacier thickness [m ice], (nbins) + running_bin_delta_thick_monthly = np.cumsum(bin_delta_thick_monthly, axis=-1) + bin_thick_monthly = running_bin_delta_thick_monthly + bin_thick_initial[:, np.newaxis] + + # --- Step 5: rebin monthly thickness --- # get surface height at the specified reference year - ref_surface_h = ds_spn.bed_h.values + ds_spn.thickness_m.sel(time=gdir.elev_change_1d['ref_dem_year']).values - - # aggregate model bin thicknesses as desired + ref_surface_height = ds_spn.bed_h.values + ds_spn.thickness_m.sel(time=gdir.elev_change_1d['ref_dem_year']).values + # aggregate model bin thicknesses with warnings.catch_warnings(): warnings.filterwarnings('ignore') - h_monthly = np.column_stack( + bin_thick_monthly = np.column_stack( [ binned_statistic( - x=ref_surface_h, values=x, statistic=np.nanmean, bins=gdir.elev_change_1d['bin_edges'] + x=ref_surface_height, + values=x, + statistic=np.nanmean, + bins=gdir.elev_change_1d['bin_edges'], )[0] - for x in h_monthly.T + for x in bin_thick_monthly.T ] ) - # interpolate over any empty bins - h_monthly_ = np.column_stack([interp1d_fill_gaps(x.copy()) for x in h_monthly.T]) + bin_thick_monthly = np.column_stack([interp1d_fill_gaps(x.copy()) for x in bin_thick_monthly.T]) - # difference desired time steps, return np.nan array for any missing inds - dh = np.column_stack( + # --- Step 6: calculate elevation change --- + elev_change_1d = np.column_stack( [ - h_monthly_[:, j] - h_monthly_[:, i] - if (i is not None and j is not None and 0 <= i < h_monthly_.shape[1] and 0 <= j < h_monthly_.shape[1]) - else np.full(h_monthly_.shape[0], np.nan) - for i, j in gdir.elev_change_1d['model2obs_inds_map'] + bin_thick_monthly[:, tup[1]] - bin_thick_monthly[:, tup[0]] + if tup[0] is not None and tup[1] is not None + else np.full(bin_thick_monthly.shape[0], np.nan) + for tup in gdir.elev_change_1d['model2obs_inds_map'] ] ) - return dh, ds_spn.dis_along_flowline.values, area + + return elev_change_1d, ds_spn.dis_along_flowline.values, area def loss_with_penalty(x, obs, mod, threshold=100, weight=1.0): @@ -286,7 +299,7 @@ def _objective(**kwargs): for start, end in gd.elev_change_1d['dates'] ] - dh_hat, dist, bin_area = get_elev_change_1d_hat(gd) + dh_hat, dist, bin_area = calculate_elev_change_1d(gd) dhdt_hat = dh_hat / gd.elev_change_1d['nyrs'] # plot binned surface area diff --git a/pygem/class_climate.py b/pygem/class_climate.py index 5cb81865..de069ee4 100755 --- a/pygem/class_climate.py +++ b/pygem/class_climate.py @@ -166,6 +166,12 @@ def __init__(self, name=str(), sim_climate_scenario=str(), realization=None): # Set parameters for ERA5, ERA-Interim, and CMIP5 netcdf files if self.name == 'ERA5': + # Ensure if using daily data, then including leap years + if pygem_prms['time']['timestep'] == 'daily': + assert pygem_prms['time']['option_leapyear'] == 1, ( + 'option_leapyear must be set to 1 if using daily ERA5 data' + ) + # Variable names self.temp_vn = 't2m' self.tempstd_vn = 't2m_std' @@ -281,6 +287,11 @@ def importGCMfxnearestneighbor_xarray(self, filename, vn, main_glac_rgi): # Import netcdf file data = xr.open_dataset(self.fx_fp + filename) glac_variable = np.zeros(main_glac_rgi.shape[0]) + + # convert longitude from -180—180 to 0—360 + if data.longitude.min() < 0: + data = data.assign_coords(longitude=(data.longitude % 360)) + # If time dimension included, then set the time index (required for ERA Interim, but not for CMIP5 or COAWST) if 'time' in data[vn].coords: time_idx = 0 @@ -599,8 +610,8 @@ def importGCMvarnearestneighbor_xarray( print('Check units of precipitation from GCM is meters per day.') if self.timestep == 'monthly' and self.name != 'COAWST': # Convert from meters per day to meters per month (COAWST data already 'monthly accumulated precipitation') - if 'daysinmonth' in dates_table.columns: - glac_variable_series = glac_variable_series * dates_table['daysinmonth'].values[np.newaxis, :] + if 'days_in_step' in dates_table.columns: + glac_variable_series = glac_variable_series * dates_table['days_in_step'].values[np.newaxis, :] elif vn != self.lr_vn: print('Check units of air temperature or precipitation') diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index b4ca6bdb..16f2e910 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -48,7 +48,6 @@ def __init__( inplace=False, debug=True, option_areaconstant=False, - spinupyears=0, constantarea_years=0, **kwargs, ): @@ -87,12 +86,10 @@ def __init__( ) self.option_areaconstant = option_areaconstant self.constantarea_years = constantarea_years - self.spinupyears = spinupyears self.glac_idx_initial = [fl.thick.nonzero()[0] for fl in flowlines] - self.y0 = 0 + self.y0 = y0 self.is_tidewater = is_tidewater self.water_level = water_level - # widths_t0 = flowlines[0].widths_m # area_v1 = widths_t0 * flowlines[0].dx_meter # print('area v1:', area_v1.sum()) @@ -326,8 +323,13 @@ def run_until_and_store(self, y1, run_path=None, diag_path=None, store_monthly_s def updategeometry(self, year, debug=False): """Update geometry for a given year""" + # get year index + year_idx = self.mb_model.get_year_index(year) + # get time step indices - note indexing should be [t_start:t_stop+1] to include final step in year + t_start, t_stop = self.mb_model.get_step_inds(year) if debug: - print('year:', year) + print('year:', year, f'({year_idx})') + print('time steps:', f'[{t_start}, {t_stop}]') # Loop over flowlines for fl_id, fl in enumerate(self.fls): @@ -338,8 +340,8 @@ def updategeometry(self, year, debug=False): width_t0 = self.fls[fl_id].widths_m.copy() # CONSTANT AREAS - # Mass redistribution ignored for calibration and spinup years (glacier properties constant) - if (self.option_areaconstant) or (year < self.spinupyears) or (year < self.constantarea_years): + # Mass redistribution ignored for calibration years (glacier properties constant) + if (self.option_areaconstant) or (year < self.y0 + self.constantarea_years): # run mass balance glac_bin_massbalclim_annual = self.mb_model.get_annual_mb( heights, fls=self.fls, fl_id=fl_id, year=year, debug=False @@ -380,7 +382,7 @@ def updategeometry(self, year, debug=False): # If frontal ablation more than bin volume, remove entire bin if fa_m3 > vol_last: # Record frontal ablation (m3 w.e.) in mass balance model for output - self.mb_model.glac_bin_frontalablation[last_idx, int(12 * (year + 1) - 1)] = ( + self.mb_model.glac_bin_frontalablation[last_idx, t_stop + 1] = ( vol_last * pygem_prms['constants']['density_ice'] / pygem_prms['constants']['density_water'] @@ -397,7 +399,7 @@ def updategeometry(self, year, debug=False): section_t0[last_idx] = section_t0[last_idx] - fa_m3 / fl.dx_meter self.fls[fl_id].section = section_t0 # Record frontal ablation(m3 w.e.) - self.mb_model.glac_bin_frontalablation[last_idx, int(12 * (year + 1) - 1)] = ( + self.mb_model.glac_bin_frontalablation[last_idx, t_stop + 1] = ( fa_m3 * pygem_prms['constants']['density_ice'] / pygem_prms['constants']['density_water'] @@ -433,11 +435,8 @@ def updategeometry(self, year, debug=False): heights, fls=self.fls, fl_id=fl_id, year=year, debug=False ) sec_in_year = ( - self.mb_model.dates_table.loc[12 * year : 12 * (year + 1) - 1, 'daysinmonth'].values.sum() - * 24 - * 3600 + self.mb_model.dates_table.iloc[t_start : t_stop + 1]['days_in_step'].values.sum() * 24 * 3600 ) - # print(' volume change [m3]:', (glac_bin_massbalclim_annual * sec_in_year * # (width_t0 * fl.dx_meter)).sum()) # print(glac_bin_masssbalclim_annual) @@ -469,14 +468,13 @@ def updategeometry(self, year, debug=False): # Record glacier properties (volume [m3], area [m2], thickness [m], width [km]) # record the next year's properties as well # 'year + 1' used so the glacier properties are consistent with mass balance computations - year = int(year) # required to ensure proper indexing with run_until_and_store (10/21/2020) glacier_area = fl.widths_m * fl.dx_meter glacier_area[fl.thick == 0] = 0 - self.mb_model.glac_bin_area_annual[:, year + 1] = glacier_area - self.mb_model.glac_bin_icethickness_annual[:, year + 1] = fl.thick - self.mb_model.glac_bin_width_annual[:, year + 1] = fl.widths_m - self.mb_model.glac_wide_area_annual[year + 1] = glacier_area.sum() - self.mb_model.glac_wide_volume_annual[year + 1] = (fl.section * fl.dx_meter).sum() + self.mb_model.glac_bin_area_annual[:, year_idx + 1] = glacier_area + self.mb_model.glac_bin_icethickness_annual[:, year_idx + 1] = fl.thick + self.mb_model.glac_bin_width_annual[:, year_idx + 1] = fl.widths_m + self.mb_model.glac_wide_area_annual[year_idx + 1] = glacier_area.sum() + self.mb_model.glac_wide_volume_annual[year_idx + 1] = (fl.section * fl.dx_meter).sum() # %% ----- FRONTAL ABLATION ----- def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, calving_k=None, debug=False): diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 60ebdb1f..20425eb7 100644 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -51,6 +51,8 @@ def __init__( Parameters ---------- + gdir : class object + Glacier directory using OGGM's structure for a single glacier modelprms : dict Model parameters dictionary (lrgcm, lrglac, precfactor, precgrad, ddfsnow, ddfice, tempsnow, tempchange) glacier_rgi_table : pd.Series @@ -63,6 +65,17 @@ def __init__( option to turn on print statements for development or debugging of code debug_refreeze : Boolean option to turn on print statements for development/debugging of refreezing code + fls : class object + flowlines using OGGM's structure for a glacier + fl_id : integer + flowline id associated with the specific flowlines to index + heights : np.array + elevation bins + inversion_filter : Boolean + option to turn on/off an inversion filter that forces the mass balance profile to be + the same or more positive with increasing elevation + ignore_debris : Boolean + option to ignore the sub-debris melt enhancement factors """ if debug: print('\n\nDEBUGGING MASS BALANCE FUNCTION\n\n') @@ -102,57 +115,57 @@ def __init__( # Variables to store (consider storing in xarray) nbins = self.glacier_area_initial.shape[0] - self.nmonths = self.glacier_gcm_temp.shape[0] + self.nsteps = self.glacier_gcm_temp.shape[0] self.years = sorted(set(self.dates_table.year.values)) self.nyears = len(self.years) # create mapper to get the appropriate year index for child functions self.year_to_index = {year: idx for idx, year in enumerate(self.years)} - self.bin_temp = np.zeros((nbins, self.nmonths)) - self.bin_prec = np.zeros((nbins, self.nmonths)) - self.bin_acc = np.zeros((nbins, self.nmonths)) - self.bin_refreezepotential = np.zeros((nbins, self.nmonths)) - self.bin_refreeze = np.zeros((nbins, self.nmonths)) - self.bin_meltglac = np.zeros((nbins, self.nmonths)) - self.bin_meltsnow = np.zeros((nbins, self.nmonths)) - self.bin_melt = np.zeros((nbins, self.nmonths)) - self.bin_snowpack = np.zeros((nbins, self.nmonths)) - self.snowpack_remaining = np.zeros((nbins, self.nmonths)) - self.glac_bin_refreeze = np.zeros((nbins, self.nmonths)) - self.glac_bin_melt = np.zeros((nbins, self.nmonths)) - self.glac_bin_frontalablation = np.zeros((nbins, self.nmonths)) - self.glac_bin_snowpack = np.zeros((nbins, self.nmonths)) - self.glac_bin_massbalclim = np.zeros((nbins, self.nmonths)) + self.bin_temp = np.zeros((nbins, self.nsteps)) + self.bin_prec = np.zeros((nbins, self.nsteps)) + self.bin_acc = np.zeros((nbins, self.nsteps)) + self.bin_refreezepotential = np.zeros((nbins, self.nsteps)) + self.bin_refreeze = np.zeros((nbins, self.nsteps)) + self.bin_meltglac = np.zeros((nbins, self.nsteps)) + self.bin_meltsnow = np.zeros((nbins, self.nsteps)) + self.bin_melt = np.zeros((nbins, self.nsteps)) + self.bin_snowpack = np.zeros((nbins, self.nsteps)) + self.snowpack_remaining = np.zeros((nbins, self.nsteps)) + self.glac_bin_refreeze = np.zeros((nbins, self.nsteps)) + self.glac_bin_melt = np.zeros((nbins, self.nsteps)) + self.glac_bin_frontalablation = np.zeros((nbins, self.nsteps)) + self.glac_bin_snowpack = np.zeros((nbins, self.nsteps)) + self.glac_bin_massbalclim = np.zeros((nbins, self.nsteps)) self.glac_bin_massbalclim_annual = np.zeros((nbins, self.nyears)) self.glac_bin_surfacetype_annual = np.zeros((nbins, self.nyears + 1)) self.glac_bin_area_annual = np.zeros((nbins, self.nyears + 1)) self.glac_bin_icethickness_annual = np.zeros((nbins, self.nyears + 1)) # Needed for MassRedistributionCurves self.glac_bin_width_annual = np.zeros((nbins, self.nyears + 1)) # Needed for MassRedistributionCurves - self.offglac_bin_prec = np.zeros((nbins, self.nmonths)) - self.offglac_bin_melt = np.zeros((nbins, self.nmonths)) - self.offglac_bin_refreeze = np.zeros((nbins, self.nmonths)) - self.offglac_bin_snowpack = np.zeros((nbins, self.nmonths)) + self.offglac_bin_prec = np.zeros((nbins, self.nsteps)) + self.offglac_bin_melt = np.zeros((nbins, self.nsteps)) + self.offglac_bin_refreeze = np.zeros((nbins, self.nsteps)) + self.offglac_bin_snowpack = np.zeros((nbins, self.nsteps)) self.offglac_bin_area_annual = np.zeros((nbins, self.nyears + 1)) - self.glac_wide_temp = np.zeros(self.nmonths) - self.glac_wide_prec = np.zeros(self.nmonths) - self.glac_wide_acc = np.zeros(self.nmonths) - self.glac_wide_refreeze = np.zeros(self.nmonths) - self.glac_wide_melt = np.zeros(self.nmonths) - self.glac_wide_frontalablation = np.zeros(self.nmonths) - self.glac_wide_massbaltotal = np.zeros(self.nmonths) - self.glac_wide_runoff = np.zeros(self.nmonths) - self.glac_wide_snowline = np.zeros(self.nmonths) + self.glac_wide_temp = np.zeros(self.nsteps) + self.glac_wide_prec = np.zeros(self.nsteps) + self.glac_wide_acc = np.zeros(self.nsteps) + self.glac_wide_refreeze = np.zeros(self.nsteps) + self.glac_wide_melt = np.zeros(self.nsteps) + self.glac_wide_frontalablation = np.zeros(self.nsteps) + self.glac_wide_massbaltotal = np.zeros(self.nsteps) + self.glac_wide_runoff = np.zeros(self.nsteps) + self.glac_wide_snowline = np.zeros(self.nsteps) self.glac_wide_area_annual = np.zeros(self.nyears + 1) self.glac_wide_volume_annual = np.zeros(self.nyears + 1) self.glac_wide_volume_change_ignored_annual = np.zeros(self.nyears) self.glac_wide_ELA_annual = np.zeros(self.nyears + 1) - self.offglac_wide_prec = np.zeros(self.nmonths) - self.offglac_wide_refreeze = np.zeros(self.nmonths) - self.offglac_wide_melt = np.zeros(self.nmonths) - self.offglac_wide_snowpack = np.zeros(self.nmonths) - self.offglac_wide_runoff = np.zeros(self.nmonths) + self.offglac_wide_prec = np.zeros(self.nsteps) + self.offglac_wide_refreeze = np.zeros(self.nsteps) + self.offglac_wide_melt = np.zeros(self.nsteps) + self.offglac_wide_snowpack = np.zeros(self.nsteps) + self.offglac_wide_runoff = np.zeros(self.nsteps) - self.dayspermonth = self.dates_table['daysinmonth'].values + self.days_in_step = self.dates_table['days_in_step'].values self.surfacetype_ddf = np.zeros((nbins)) # Surface type DDF dictionary (manipulate this function for calibration or for each glacier) @@ -182,13 +195,9 @@ def __init__( # refrezee cold content or "potential" refreeze self.rf_cold = np.zeros(nbins) # layer temp of each elev bin for present time step - self.te_rf = np.zeros((pygem_prms['mb']['HH2015_rf_opts']['rf_layers'], nbins, self.nmonths)) + self.te_rf = np.zeros((pygem_prms['mb']['HH2015_rf_opts']['rf_layers'], nbins, self.nsteps)) # layer temp of each elev bin for previous time step - self.tl_rf = np.zeros((pygem_prms['mb']['HH2015_rf_opts']['rf_layers'], nbins, self.nmonths)) - - # Sea level for marine-terminating glaciers - self.sea_level = 0 - rgi_region = int(glacier_rgi_table.RGIId.split('-')[1].split('.')[0]) + self.tl_rf = np.zeros((pygem_prms['mb']['HH2015_rf_opts']['rf_layers'], nbins, self.nsteps)) def get_year_index(self, year): """ @@ -198,6 +207,18 @@ def get_year_index(self, year): assert year in self.years, f'{year} not found in model dates table' return self.year_to_index[year] + def get_step_inds(self, year): + """ + Get time step start and stop indices for the specified model year. + + Both indices are inclusive, so standard Python slicing would use + [t_start : t_stop + 1]. + """ + step_idxs = np.where(self.dates_table.year == int(year))[0] + if step_idxs.size == 0: + raise ValueError(f'Year {year} not found in dates_table.') + return step_idxs[0], step_idxs[-1] + def get_annual_mb( self, heights, @@ -223,13 +244,16 @@ def get_annual_mb( mb : np.array mass balance for each bin [m ice per second] """ + + # assertion to only run with calendar years + assert pygem_prms['climate']['sim_wateryear'] == 'calendar', ( + 'This function is not set up yet to handle non-calendar years' + ) + # get year index year_idx = self.get_year_index(year) - # get start step for 0th month of specified year - year_start_month_idx = 12 * year_idx - # get stop step for specified year - # note, this is 1 greater than the final month which to include - python indexing will not include this month, final month to include of given year is <12*(year_idx+1)-1> - year_stop_month_idx = 12 * (year_idx + 1) + # get time step indices - note indexing should be [t_start:t_stop+1] to include final step in year + t_start, t_stop = self.get_step_inds(year) fl = fls[fl_id] @@ -256,10 +280,10 @@ def get_annual_mb( glac_idx_t0 = glacier_area_t0.nonzero()[0] nbins = heights.shape[0] - nmonths = self.glacier_gcm_temp.shape[0] + nsteps = self.glacier_gcm_temp.shape[0] # Local variables - bin_precsnow = np.zeros((nbins, nmonths)) + bin_precsnow = np.zeros((nbins, nsteps)) # Refreezing specific layers if pygem_prms['mb']['option_refreezing'] == 'HH2015' and year_idx == 0: @@ -269,8 +293,6 @@ def get_annual_mb( refreeze_potential = np.zeros(nbins) if self.glacier_area_initial.sum() > 0: - # if len(glac_idx_t0) > 0: - # Surface type [0=off-glacier, 1=ice, 2=snow, 3=firn, 4=debris] if year_idx == 0: self.surfacetype, self.firnline_idx = self._surfacetypebinsinitial(self.heights) @@ -281,507 +303,482 @@ def get_annual_mb( self.offglac_bin_area_annual[:, year_idx] = glacier_area_initial - glacier_area_t0 offglac_idx = np.where(self.offglac_bin_area_annual[:, year_idx] > 0)[0] - # Functions currently set up for monthly timestep - # only compute mass balance while glacier exists - if pygem_prms['time']['timestep'] == 'monthly': - # if (pygem_prms['time']['timestep'] == 'monthly') and (glac_idx_t0.shape[0] != 0): - - # AIR TEMPERATURE: Downscale the gcm temperature [deg C] to each bin - # Downscale using gcm and glacier lapse rates - # T_bin = T_gcm + lr_gcm * (z_ref - z_gcm) + lr_glac * (z_bin - z_ref) + tempchange - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] = ( - self.glacier_gcm_temp[year_start_month_idx:year_stop_month_idx] - + self.glacier_gcm_lrgcm[year_start_month_idx:year_stop_month_idx] - * ( - self.glacier_rgi_table.loc[pygem_prms['mb']['option_elev_ref_downscale']] - - self.glacier_gcm_elev + # AIR TEMPERATURE: Downscale the gcm temperature [deg C] to each bin + # Downscale using gcm and glacier lapse rates + # T_bin = T_gcm + lr_gcm * (z_ref - z_gcm) + lr_glac * (z_bin - z_ref) + tempchange + self.bin_temp[:, t_start : t_stop + 1] = ( + self.glacier_gcm_temp[t_start : t_stop + 1] + + self.glacier_gcm_lrgcm[t_start : t_stop + 1] + * (self.glacier_rgi_table.loc[pygem_prms['mb']['option_elev_ref_downscale']] - self.glacier_gcm_elev) + + self.glacier_gcm_lrglac[t_start : t_stop + 1] + * (heights - self.glacier_rgi_table.loc[pygem_prms['mb']['option_elev_ref_downscale']])[:, np.newaxis] + + self.modelprms['tbias'] + ) + + # PRECIPITATION/ACCUMULATION: Downscale the precipitation (liquid and solid) to each bin + # Precipitation using precipitation factor and precipitation gradient + # P_bin = P_gcm * prec_factor * (1 + prec_grad * (z_bin - z_ref)) + bin_precsnow[:, t_start : t_stop + 1] = ( + self.glacier_gcm_prec[t_start : t_stop + 1] + * self.modelprms['kp'] + * ( + 1 + + self.modelprms['precgrad'] + * (heights - self.glacier_rgi_table.loc[pygem_prms['mb']['option_elev_ref_downscale']]) + )[:, np.newaxis] + ) + # Option to adjust prec of uppermost 25% of glacier for wind erosion and reduced moisture content + if pygem_prms['mb']['option_preclimit'] == 1: + # Elevation range based on all flowlines + raw_min_elev = [] + raw_max_elev = [] + if len(fl.surface_h[fl.widths_m > 0]): + raw_min_elev.append(fl.surface_h[fl.widths_m > 0].min()) + raw_max_elev.append(fl.surface_h[fl.widths_m > 0].max()) + elev_range = np.max(raw_max_elev) - np.min(raw_min_elev) + elev_75 = np.min(raw_min_elev) + 0.75 * (elev_range) + + # If elevation range > 1000 m, apply corrections to uppermost 25% of glacier (Huss and Hock, 2015) + if elev_range > 1000: + # Indices of upper 25% + glac_idx_upper25 = glac_idx_t0[heights[glac_idx_t0] >= elev_75] + # Exponential decay according to elevation difference from the 75% elevation + # prec_upper25 = prec * exp(-(elev_i - elev_75%)/(elev_max- - elev_75%)) + # height at 75% of the elevation + height_75 = heights[glac_idx_upper25].min() + glac_idx_75 = np.where(heights == height_75)[0][0] + # exponential decay + bin_precsnow[glac_idx_upper25, t_start : t_stop + 1] = ( + bin_precsnow[glac_idx_75, t_start : t_stop + 1] + * np.exp( + -1 + * (heights[glac_idx_upper25] - height_75) + / (heights[glac_idx_upper25].max() - heights[glac_idx_upper25].min()) + )[:, np.newaxis] + ) + # Precipitation cannot be less than 87.5% of the maximum accumulation elsewhere on the glacier + # compute max values for each step over glac_idx_t0, compare, and replace if needed + max_values = np.tile( + 0.875 * np.max(bin_precsnow[glac_idx_t0, t_start : t_stop + 1], axis=0), + (8, 1), ) - + self.glacier_gcm_lrglac[year_start_month_idx:year_stop_month_idx] - * (heights - self.glacier_rgi_table.loc[pygem_prms['mb']['option_elev_ref_downscale']])[ - :, np.newaxis + uncorrected_values = bin_precsnow[glac_idx_upper25, t_start : t_stop + 1] + corrected_values = np.max(np.stack([uncorrected_values, max_values], axis=0), axis=0) + bin_precsnow[glac_idx_upper25, t_start : t_stop + 1] = corrected_values + + # Separate total precipitation into liquid (bin_prec) and solid (bin_acc) + if pygem_prms['mb']['option_accumulation'] == 1: + # if temperature above threshold, then rain + ( + self.bin_prec[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] ] - + self.modelprms['tbias'] - ) - - # PRECIPITATION/ACCUMULATION: Downscale the precipitation (liquid and solid) to each bin - # Precipitation using precipitation factor and precipitation gradient - # P_bin = P_gcm * prec_factor * (1 + prec_grad * (z_bin - z_ref)) - bin_precsnow[:, year_start_month_idx:year_stop_month_idx] = ( - self.glacier_gcm_prec[year_start_month_idx:year_stop_month_idx] - * self.modelprms['kp'] - * ( - 1 - + self.modelprms['precgrad'] - * (heights - self.glacier_rgi_table.loc[pygem_prms['mb']['option_elev_ref_downscale']]) - )[:, np.newaxis] + ) = bin_precsnow[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] + ] + # if temperature below threshold, then snow + ( + self.bin_acc[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] + ] + ) = bin_precsnow[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] + ] + elif pygem_prms['mb']['option_accumulation'] == 2: + # if temperature between min/max, then mix of snow/rain using linear relationship between min/max + self.bin_prec[:, t_start : t_stop + 1] = ( + 0.5 + (self.bin_temp[:, t_start : t_stop + 1] - self.modelprms['tsnow_threshold']) / 2 + ) * bin_precsnow[:, t_start : t_stop + 1] + self.bin_acc[:, t_start : t_stop + 1] = ( + bin_precsnow[:, t_start : t_stop + 1] - self.bin_prec[:, t_start : t_stop + 1] ) - # Option to adjust prec of uppermost 25% of glacier for wind erosion and reduced moisture content - if pygem_prms['mb']['option_preclimit'] == 1: - # Elevation range based on all flowlines - raw_min_elev = [] - raw_max_elev = [] - if len(fl.surface_h[fl.widths_m > 0]): - raw_min_elev.append(fl.surface_h[fl.widths_m > 0].min()) - raw_max_elev.append(fl.surface_h[fl.widths_m > 0].max()) - elev_range = np.max(raw_max_elev) - np.min(raw_min_elev) - elev_75 = np.min(raw_min_elev) + 0.75 * (elev_range) - - # If elevation range > 1000 m, apply corrections to uppermost 25% of glacier (Huss and Hock, 2015) - if elev_range > 1000: - # Indices of upper 25% - glac_idx_upper25 = glac_idx_t0[heights[glac_idx_t0] >= elev_75] - # Exponential decay according to elevation difference from the 75% elevation - # prec_upper25 = prec * exp(-(elev_i - elev_75%)/(elev_max- - elev_75%)) - # height at 75% of the elevation - height_75 = heights[glac_idx_upper25].min() - glac_idx_75 = np.where(heights == height_75)[0][0] - # exponential decay - bin_precsnow[glac_idx_upper25, year_start_month_idx:year_stop_month_idx] = ( - bin_precsnow[glac_idx_75, year_start_month_idx:year_stop_month_idx] - * np.exp( - -1 - * (heights[glac_idx_upper25] - height_75) - / (heights[glac_idx_upper25].max() - heights[glac_idx_upper25].min()) - )[:, np.newaxis] - ) - # Precipitation cannot be less than 87.5% of the maximum accumulation elsewhere on the glacier - for month in range(0, 12): - bin_precsnow[ - glac_idx_upper25[ - ( - bin_precsnow[glac_idx_upper25, month] - < 0.875 * bin_precsnow[glac_idx_t0, month].max() - ) - & (bin_precsnow[glac_idx_upper25, month] != 0) - ], - month, - ] = 0.875 * bin_precsnow[glac_idx_t0, month].max() - - # Separate total precipitation into liquid (bin_prec) and solid (bin_acc) - if pygem_prms['mb']['option_accumulation'] == 1: - # if temperature above threshold, then rain - ( - self.bin_prec[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] - > self.modelprms['tsnow_threshold'] - ] - ) = bin_precsnow[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] > self.modelprms['tsnow_threshold'] + # if temperature above maximum threshold, then all rain + ( + self.bin_prec[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] + 1 ] - # if temperature below threshold, then snow - ( - self.bin_acc[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] - <= self.modelprms['tsnow_threshold'] - ] - ) = bin_precsnow[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] <= self.modelprms['tsnow_threshold'] + ) = bin_precsnow[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] + 1 + ] + ( + self.bin_acc[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] + 1 ] - elif pygem_prms['mb']['option_accumulation'] == 2: - # if temperature between min/max, then mix of snow/rain using linear relationship between min/max - self.bin_prec[:, year_start_month_idx:year_stop_month_idx] = ( - 0.5 - + ( - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] - - self.modelprms['tsnow_threshold'] - ) - / 2 - ) * bin_precsnow[:, year_start_month_idx:year_stop_month_idx] - self.bin_acc[:, year_start_month_idx:year_stop_month_idx] = ( - bin_precsnow[:, year_start_month_idx:year_stop_month_idx] - - self.bin_prec[:, year_start_month_idx:year_stop_month_idx] - ) - # if temperature above maximum threshold, then all rain - ( - self.bin_prec[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] - > self.modelprms['tsnow_threshold'] + 1 - ] - ) = bin_precsnow[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] - > self.modelprms['tsnow_threshold'] + 1 + ) = 0 + # if temperature below minimum threshold, then all snow + ( + self.bin_acc[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] - 1 ] - ( - self.bin_acc[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] - > self.modelprms['tsnow_threshold'] + 1 - ] - ) = 0 - # if temperature below minimum threshold, then all snow - ( - self.bin_acc[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] - <= self.modelprms['tsnow_threshold'] - 1 - ] - ) = bin_precsnow[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] - <= self.modelprms['tsnow_threshold'] - 1 + ) = bin_precsnow[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] - 1 + ] + ( + self.bin_prec[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] - 1 ] - ( - self.bin_prec[:, year_start_month_idx:year_stop_month_idx][ - self.bin_temp[:, year_start_month_idx:year_stop_month_idx] - <= self.modelprms['tsnow_threshold'] - 1 - ] - ) = 0 - - # ENTER MONTHLY LOOP (monthly loop required since surface type changes) - for month in range(0, 12): - # Step is the position as a function of year and month, which improves readability - step = year_start_month_idx + month - - # ACCUMULATION, MELT, REFREEZE, AND CLIMATIC MASS BALANCE - # Snowpack [m w.e.] = snow remaining + new snow - if step == 0: - self.bin_snowpack[:, step] = self.bin_acc[:, step] - else: - self.bin_snowpack[:, step] = self.snowpack_remaining[:, step - 1] + self.bin_acc[:, step] - - # MELT [m w.e.] - # energy available for melt [degC day] - if pygem_prms['mb']['option_ablation'] == 1: - # option 1: energy based on monthly temperature - melt_energy_available = self.bin_temp[:, step] * self.dayspermonth[step] - melt_energy_available[melt_energy_available < 0] = 0 - elif pygem_prms['mb']['option_ablation'] == 2: - # Seed randomness for repeatability, but base it on step to ensure the daily variability is not - # the same for every single time step - np.random.seed(step) - # option 2: monthly temperature superimposed with daily temperature variability - # daily temperature variation in each bin for the monthly timestep - bin_tempstd_daily = np.repeat( - np.random.normal( - loc=0, - scale=self.glacier_gcm_tempstd[step], - size=self.dayspermonth[step], - ).reshape(1, self.dayspermonth[step]), - heights.shape[0], - axis=0, - ) - # daily temperature in each bin for the monthly timestep - bin_temp_daily = self.bin_temp[:, step][:, np.newaxis] + bin_tempstd_daily - # remove negative values - bin_temp_daily[bin_temp_daily < 0] = 0 - # Energy available for melt [degC day] = sum of daily energy available - melt_energy_available = bin_temp_daily.sum(axis=1) - # SNOW MELT [m w.e.] - self.bin_meltsnow[:, step] = self.surfacetype_ddf_dict[2] * melt_energy_available - # snow melt cannot exceed the snow depth - self.bin_meltsnow[self.bin_meltsnow[:, step] > self.bin_snowpack[:, step], step] = ( - self.bin_snowpack[self.bin_meltsnow[:, step] > self.bin_snowpack[:, step], step] + ) = 0 + + # ENTER TIME STEP LOOP (loop is required since surface type changes) + for step in range(t_start, t_stop): + # ACCUMULATION, MELT, REFREEZE, AND CLIMATIC MASS BALANCE + # Snowpack [m w.e.] = snow remaining + new snow + if step == 0: + self.bin_snowpack[:, step] = self.bin_acc[:, step] + else: + self.bin_snowpack[:, step] = self.snowpack_remaining[:, step - 1] + self.bin_acc[:, step] + + # MELT [m w.e.] + # energy available for melt [degC day] + if pygem_prms['mb']['option_ablation'] == 1: + # option 1: energy based on temperature + melt_energy_available = self.bin_temp[:, step] * self.days_in_step[step] + # assert 1==0, 'here is where we need to change to days per step' + melt_energy_available[melt_energy_available < 0] = 0 + + elif pygem_prms['mb']['option_ablation'] == 2: + assert pygem_prms['time']['timestep'] != 'daily', ( + 'Option 2 for ablation should not be used with daily data' ) - # GLACIER MELT (ice and firn) [m w.e.] - # energy remaining after snow melt [degC day] - melt_energy_available = ( - melt_energy_available - self.bin_meltsnow[:, step] / self.surfacetype_ddf_dict[2] + + # option 2: monthly temperature superimposed with daily temperature variability + # daily temperature variation in each bin for the monthly timestep + # Seed randomness for repeatability, but base it on step to ensure the daily variability is not + # the same for every single time step + np.random.seed(step) + + bin_tempstd_daily = np.repeat( + np.random.normal( + loc=0, + scale=self.glacier_gcm_tempstd[step], + size=self.days_in_step[step], + ).reshape(1, self.days_in_step[step]), + heights.shape[0], + axis=0, ) - # remove low values of energy available caused by rounding errors in the step above - melt_energy_available[abs(melt_energy_available) < pygem_prms['constants']['tolerance']] = 0 - # DDF based on surface type [m w.e. degC-1 day-1] - for surfacetype_idx in self.surfacetype_ddf_dict: - self.surfacetype_ddf[self.surfacetype == surfacetype_idx] = self.surfacetype_ddf_dict[ - surfacetype_idx - ] - # Debris enhancement factors in ablation area (debris in accumulation area would submerge) - if surfacetype_idx == 1 and pygem_prms['mb']['include_debris']: - self.surfacetype_ddf[self.surfacetype == 1] = ( - self.surfacetype_ddf[self.surfacetype == 1] * self.debris_ed[self.surfacetype == 1] - ) - self.bin_meltglac[glac_idx_t0, step] = ( - self.surfacetype_ddf[glac_idx_t0] * melt_energy_available[glac_idx_t0] + # daily temperature in each bin for the given timestep + bin_temp_daily = self.bin_temp[:, step][:, np.newaxis] + bin_tempstd_daily + # remove negative values + bin_temp_daily[bin_temp_daily < 0] = 0 + # Energy available for melt [degC day] = sum of daily energy available + melt_energy_available = bin_temp_daily.sum(axis=1) + # SNOW MELT [m w.e.] + self.bin_meltsnow[:, step] = self.surfacetype_ddf_dict[2] * melt_energy_available + # snow melt cannot exceed the snow depth + self.bin_meltsnow[self.bin_meltsnow[:, step] > self.bin_snowpack[:, step], step] = self.bin_snowpack[ + self.bin_meltsnow[:, step] > self.bin_snowpack[:, step], step + ] + # GLACIER MELT (ice and firn) [m w.e.] + # energy remaining after snow melt [degC day] + melt_energy_available = ( + melt_energy_available - self.bin_meltsnow[:, step] / self.surfacetype_ddf_dict[2] + ) + # remove low values of energy available caused by rounding errors in the step above + melt_energy_available[abs(melt_energy_available) < pygem_prms['constants']['tolerance']] = 0 + # DDF based on surface type [m w.e. degC-1 day-1] + for surfacetype_idx in self.surfacetype_ddf_dict: + self.surfacetype_ddf[self.surfacetype == surfacetype_idx] = self.surfacetype_ddf_dict[ + surfacetype_idx + ] + # Debris enhancement factors in ablation area (debris in accumulation area would submerge) + if surfacetype_idx == 1 and pygem_prms['mb']['include_debris']: + self.surfacetype_ddf[self.surfacetype == 1] = ( + self.surfacetype_ddf[self.surfacetype == 1] * self.debris_ed[self.surfacetype == 1] + ) + self.bin_meltglac[glac_idx_t0, step] = ( + self.surfacetype_ddf[glac_idx_t0] * melt_energy_available[glac_idx_t0] + ) + # TOTAL MELT (snow + glacier) + # off-glacier need to include melt of refreeze because there are no glacier dynamics, + # but on-glacier do not need to account for this (simply assume refreeze has same surface type) + self.bin_melt[:, step] = self.bin_meltglac[:, step] + self.bin_meltsnow[:, step] + + # REFREEZING + if pygem_prms['mb']['option_refreezing'] == 'HH2015': + if step > 0: + self.tl_rf[:, :, step] = self.tl_rf[:, :, step - 1] + self.te_rf[:, :, step] = self.te_rf[:, :, step - 1] + + # Refreeze based on heat conduction approach (Huss and Hock 2015) + # refreeze time step (s) + rf_dt = 3600 * 24 * self.days_in_step[step] / pygem_prms['mb']['HH2015_rf_opts']['rf_dsc'] + + if pygem_prms['mb']['HH2015_rf_opts']['option_rf_limit_meltsnow'] == 1: + bin_meltlimit = self.bin_meltsnow.copy() + else: + bin_meltlimit = self.bin_melt.copy() + + # Debug lowest bin + if self.debug_refreeze: + gidx_debug = np.where(heights == heights[glac_idx_t0].min())[0] + + # Loop through each elevation bin of glacier + assert pygem_prms['time']['timestep'] != 'daily', ( + 'must remove the 12 thats hard-coded here - did not do this for HH2015 given the issue' ) - # TOTAL MELT (snow + glacier) - # off-glacier need to include melt of refreeze because there are no glacier dynamics, - # but on-glacier do not need to account for this (simply assume refreeze has same surface type) - self.bin_melt[:, step] = self.bin_meltglac[:, step] + self.bin_meltsnow[:, step] - - # REFREEZING - if pygem_prms['mb']['option_refreezing'] == 'HH2015': - if step > 0: - self.tl_rf[:, :, step] = self.tl_rf[:, :, step - 1] - self.te_rf[:, :, step] = self.te_rf[:, :, step - 1] - - # Refreeze based on heat conduction approach (Huss and Hock 2015) - # refreeze time step (s) - rf_dt = 3600 * 24 * self.dayspermonth[step] / pygem_prms['mb']['HH2015_rf_opts']['rf_dsc'] - - if pygem_prms['mb']['HH2015_rf_opts']['option_rf_limit_meltsnow'] == 1: - bin_meltlimit = self.bin_meltsnow.copy() - else: - bin_meltlimit = self.bin_melt.copy() - # Debug lowest bin - if self.debug_refreeze: - gidx_debug = np.where(heights == heights[glac_idx_t0].min())[0] + for nbin, gidx in enumerate(glac_idx_t0): + # COMPUTE HEAT CONDUCTION - BUILD COLD RESERVOIR + # If no melt, then build up cold reservoir (compute heat conduction) + if self.bin_melt[gidx, step] < pygem_prms['mb']['HH2015_rf_opts']['rf_meltcrit']: + if self.debug_refreeze and gidx == gidx_debug and step < 12: + print( + '\nMonth ' + str(self.dates_table.loc[step, 'month']), + 'Computing heat conduction', + ) - # Loop through each elevation bin of glacier - for nbin, gidx in enumerate(glac_idx_t0): - # COMPUTE HEAT CONDUCTION - BUILD COLD RESERVOIR - # If no melt, then build up cold reservoir (compute heat conduction) - if self.bin_melt[gidx, step] < pygem_prms['mb']['HH2015_rf_opts']['rf_meltcrit']: - if self.debug_refreeze and gidx == gidx_debug and step < 12: - print( - '\nMonth ' + str(self.dates_table.loc[step, 'month']), - 'Computing heat conduction', + # Set refreeze equal to 0 + self.refr[gidx] = 0 + # Loop through multiple iterations to converge on a solution + # -> this will loop through 0, 1, 2 + for h in np.arange(0, pygem_prms['mb']['HH2015_rf_opts']['rf_dsc']): + # Compute heat conduction in layers (loop through rows) + # go from 1 to rf_layers-1 to avoid indexing errors with "j-1" and "j+1" + # "j+1" is set to zero, which is fine for temperate glaciers but inaccurate for + # cold/polythermal glaciers + for j in np.arange( + 1, + pygem_prms['mb']['HH2015_rf_opts']['rf_layers'] - 1, + ): + # Assume temperature of first layer equals air temperature + # assumption probably wrong, but might still work at annual average + # Since next line uses tl_rf for all calculations, set tl_rf[0] to present mean + # monthly air temperature to ensure the present calculations are done with the + # present time step's air temperature + self.tl_rf[0, gidx, step] = self.bin_temp[gidx, step] + # Temperature for each layer + self.te_rf[j, gidx, step] = self.tl_rf[j, gidx, step] + rf_dt * self.rf_layers_k[ + j + ] / self.rf_layers_ch[j] / pygem_prms['mb']['HH2015_rf_opts'][ + 'rf_dz' + ] ** 2 * 0.5 * ( + (self.tl_rf[j - 1, gidx, step] - self.tl_rf[j, gidx, step]) + - (self.tl_rf[j, gidx, step] - self.tl_rf[j + 1, gidx, step]) ) + # Update previous time step + self.tl_rf[:, gidx, step] = self.te_rf[:, gidx, step] - # Set refreeze equal to 0 - self.refr[gidx] = 0 - # Loop through multiple iterations to converge on a solution - # -> this will loop through 0, 1, 2 - for h in np.arange(0, pygem_prms['mb']['HH2015_rf_opts']['rf_dsc']): - # Compute heat conduction in layers (loop through rows) - # go from 1 to rf_layers-1 to avoid indexing errors with "j-1" and "j+1" - # "j+1" is set to zero, which is fine for temperate glaciers but inaccurate for - # cold/polythermal glaciers - for j in np.arange( - 1, - pygem_prms['mb']['HH2015_rf_opts']['rf_layers'] - 1, - ): - # Assume temperature of first layer equals air temperature - # assumption probably wrong, but might still work at annual average - # Since next line uses tl_rf for all calculations, set tl_rf[0] to present mean - # monthly air temperature to ensure the present calculations are done with the - # present time step's air temperature - self.tl_rf[0, gidx, step] = self.bin_temp[gidx, step] - # Temperature for each layer - self.te_rf[j, gidx, step] = self.tl_rf[ - j, gidx, step - ] + rf_dt * self.rf_layers_k[j] / self.rf_layers_ch[j] / pygem_prms['mb'][ - 'HH2015_rf_opts' - ]['rf_dz'] ** 2 * 0.5 * ( - (self.tl_rf[j - 1, gidx, step] - self.tl_rf[j, gidx, step]) - - (self.tl_rf[j, gidx, step] - self.tl_rf[j + 1, gidx, step]) - ) - # Update previous time step - self.tl_rf[:, gidx, step] = self.te_rf[:, gidx, step] + if self.debug_refreeze and gidx == gidx_debug and step < 12: + print( + 'tl_rf:', + ['{:.2f}'.format(x) for x in self.tl_rf[:, gidx, step]], + ) - if self.debug_refreeze and gidx == gidx_debug and step < 12: - print( - 'tl_rf:', - ['{:.2f}'.format(x) for x in self.tl_rf[:, gidx, step]], - ) + # COMPUTE REFREEZING - TAP INTO "COLD RESERVOIR" or potential refreezing + else: + if self.debug_refreeze and gidx == gidx_debug and step < 12: + print( + '\nMonth ' + str(self.dates_table.loc[step, 'month']), + 'Computing refreeze', + ) - # COMPUTE REFREEZING - TAP INTO "COLD RESERVOIR" or potential refreezing + # Refreezing over firn surface + if (self.surfacetype[gidx] == 2) or (self.surfacetype[gidx] == 3): + nlayers = pygem_prms['mb']['HH2015_rf_opts']['rf_layers'] - 1 + # Refreezing over ice surface else: + # Approximate number of layers of snow on top of ice + smax = np.round( + ( + self.bin_snowpack[gidx, step] / (self.rf_layers_dens[0] / 1000) + + pygem_prms['mb']['HH2015_rf_opts']['pp'] + ) + / pygem_prms['mb']['HH2015_rf_opts']['rf_dz'], + 0, + ) + # if there is very little snow on the ground (SWE > 0.06 m for pp=0.3), + # then still set smax (layers) to 1 + if self.bin_snowpack[gidx, step] > 0 and smax == 0: + smax = 1 + # if no snow on the ground, then set to rf_cold to NoData value + if smax == 0: + self.rf_cold[gidx] = 0 + # if smax greater than the number of layers, set to max number of layers minus 1 + if smax > pygem_prms['mb']['HH2015_rf_opts']['rf_layers'] - 1: + smax = pygem_prms['mb']['HH2015_rf_opts']['rf_layers'] - 1 + nlayers = int(smax) + # Compute potential refreeze, "cold reservoir", from temperature in each layer + # only calculate potential refreezing first time it starts melting each year + if self.rf_cold[gidx] == 0 and self.tl_rf[:, gidx, step].min() < 0: if self.debug_refreeze and gidx == gidx_debug and step < 12: - print( - '\nMonth ' + str(self.dates_table.loc[step, 'month']), - 'Computing refreeze', + print('calculating potential refreeze from ' + str(nlayers) + ' layers') + + for j in np.arange(0, nlayers): + j += 1 + # units: (degC) * (J K-1 m-3) * (m) * (kg J-1) * (m3 kg-1) + rf_cold_layer = ( + self.tl_rf[j, gidx, step] + * self.rf_layers_ch[j] + * pygem_prms['mb']['HH2015_rf_opts']['rf_dz'] + / pygem_prms['constants']['Lh_rf'] + / pygem_prms['constants']['density_water'] ) + self.rf_cold[gidx] -= rf_cold_layer - # Refreezing over firn surface - if (self.surfacetype[gidx] == 2) or (self.surfacetype[gidx] == 3): - nlayers = pygem_prms['mb']['HH2015_rf_opts']['rf_layers'] - 1 - # Refreezing over ice surface - else: - # Approximate number of layers of snow on top of ice - smax = np.round( - ( - self.bin_snowpack[gidx, step] / (self.rf_layers_dens[0] / 1000) - + pygem_prms['mb']['HH2015_rf_opts']['pp'] - ) - / pygem_prms['mb']['HH2015_rf_opts']['rf_dz'], - 0, - ) - # if there is very little snow on the ground (SWE > 0.06 m for pp=0.3), - # then still set smax (layers) to 1 - if self.bin_snowpack[gidx, step] > 0 and smax == 0: - smax = 1 - # if no snow on the ground, then set to rf_cold to NoData value - if smax == 0: - self.rf_cold[gidx] = 0 - # if smax greater than the number of layers, set to max number of layers minus 1 - if smax > pygem_prms['mb']['HH2015_rf_opts']['rf_layers'] - 1: - smax = pygem_prms['mb']['HH2015_rf_opts']['rf_layers'] - 1 - nlayers = int(smax) - # Compute potential refreeze, "cold reservoir", from temperature in each layer - # only calculate potential refreezing first time it starts melting each year - if self.rf_cold[gidx] == 0 and self.tl_rf[:, gidx, step].min() < 0: if self.debug_refreeze and gidx == gidx_debug and step < 12: - print('calculating potential refreeze from ' + str(nlayers) + ' layers') - - for j in np.arange(0, nlayers): - j += 1 - # units: (degC) * (J K-1 m-3) * (m) * (kg J-1) * (m3 kg-1) - rf_cold_layer = ( - self.tl_rf[j, gidx, step] - * self.rf_layers_ch[j] - * pygem_prms['mb']['HH2015_rf_opts']['rf_dz'] - / pygem_prms['constants']['Lh_rf'] - / pygem_prms['constants']['density_water'] + print( + 'j:', + j, + 'tl_rf @ j:', + np.round(self.tl_rf[j, gidx, step], 2), + 'ch @ j:', + np.round(self.rf_layers_ch[j], 2), + 'rf_cold_layer @ j:', + np.round(rf_cold_layer, 2), + 'rf_cold @ j:', + np.round(self.rf_cold[gidx], 2), ) - self.rf_cold[gidx] -= rf_cold_layer - - if self.debug_refreeze and gidx == gidx_debug and step < 12: - print( - 'j:', - j, - 'tl_rf @ j:', - np.round(self.tl_rf[j, gidx, step], 2), - 'ch @ j:', - np.round(self.rf_layers_ch[j], 2), - 'rf_cold_layer @ j:', - np.round(rf_cold_layer, 2), - 'rf_cold @ j:', - np.round(self.rf_cold[gidx], 2), - ) - if self.debug_refreeze and gidx == gidx_debug and step < 12: - print('rf_cold:', np.round(self.rf_cold[gidx], 2)) - - # Compute refreezing - # If melt and liquid prec < potential refreeze, then refreeze all melt and liquid prec - if (bin_meltlimit[gidx, step] + self.bin_prec[gidx, step]) < self.rf_cold[gidx]: - self.refr[gidx] = bin_meltlimit[gidx, step] + self.bin_prec[gidx, step] - # otherwise, refreeze equals the potential refreeze - elif self.rf_cold[gidx] > 0: - self.refr[gidx] = self.rf_cold[gidx] - else: - self.refr[gidx] = 0 - - # Track the remaining potential refreeze - self.rf_cold[gidx] -= bin_meltlimit[gidx, step] + self.bin_prec[gidx, step] - # if potential refreeze consumed, set to 0 and set temperature to 0 (temperate firn) - if self.rf_cold[gidx] < 0: - self.rf_cold[gidx] = 0 - self.tl_rf[:, gidx, step] = 0 + if self.debug_refreeze and gidx == gidx_debug and step < 12: + print('rf_cold:', np.round(self.rf_cold[gidx], 2)) + + # Compute refreezing + # If melt and liquid prec < potential refreeze, then refreeze all melt and liquid prec + if (bin_meltlimit[gidx, step] + self.bin_prec[gidx, step]) < self.rf_cold[gidx]: + self.refr[gidx] = bin_meltlimit[gidx, step] + self.bin_prec[gidx, step] + # otherwise, refreeze equals the potential refreeze + elif self.rf_cold[gidx] > 0: + self.refr[gidx] = self.rf_cold[gidx] + else: + self.refr[gidx] = 0 - # Record refreeze - self.bin_refreeze[gidx, step] = self.refr[gidx] + # Track the remaining potential refreeze + self.rf_cold[gidx] -= bin_meltlimit[gidx, step] + self.bin_prec[gidx, step] + # if potential refreeze consumed, set to 0 and set temperature to 0 (temperate firn) + if self.rf_cold[gidx] < 0: + self.rf_cold[gidx] = 0 + self.tl_rf[:, gidx, step] = 0 - if self.debug_refreeze and step < 12 and gidx == gidx_debug: - print( - 'Month ' + str(self.dates_table.loc[step, 'month']), - 'Rf_cold remaining:', - np.round(self.rf_cold[gidx], 2), - 'Snow depth:', - np.round(self.bin_snowpack[glac_idx_t0[nbin], step], 2), - 'Snow melt:', - np.round(self.bin_meltsnow[glac_idx_t0[nbin], step], 2), - 'Rain:', - np.round(self.bin_prec[glac_idx_t0[nbin], step], 2), - 'Rfrz:', - np.round(self.bin_refreeze[gidx, step], 2), - ) + # Record refreeze + self.bin_refreeze[gidx, step] = self.refr[gidx] - elif pygem_prms['mb']['option_refreezing'] == 'Woodward': - # Refreeze based on annual air temperature (Woodward etal. 1997) - # R(m) = (-0.69 * Tair + 0.0096) * 1 m / 100 cm - # calculate annually and place potential refreeze in user defined month - if step % 12 == 0: - bin_temp_annual = annualweightedmean_array( - self.bin_temp[:, year_start_month_idx:year_stop_month_idx], - self.dates_table.iloc[year_start_month_idx:year_stop_month_idx, :], - ) - bin_refreezepotential_annual = (-0.69 * bin_temp_annual + 0.0096) / 100 - # Remove negative refreezing values - bin_refreezepotential_annual[bin_refreezepotential_annual < 0] = 0 - self.bin_refreezepotential[:, step] = bin_refreezepotential_annual - # Reset refreeze potential every year - if self.bin_refreezepotential[:, step].max() > 0: - refreeze_potential = self.bin_refreezepotential[:, step] - - if self.debug_refreeze: + if self.debug_refreeze and step < 12 and gidx == gidx_debug: print( - 'Year ' + str(year) + ' Month ' + str(self.dates_table.loc[step, 'month']), - 'Refreeze potential:', - np.round(refreeze_potential[glac_idx_t0[0]], 3), + 'Month ' + str(self.dates_table.loc[step, 'month']), + 'Rf_cold remaining:', + np.round(self.rf_cold[gidx], 2), 'Snow depth:', - np.round(self.bin_snowpack[glac_idx_t0[0], step], 2), + np.round(self.bin_snowpack[glac_idx_t0[nbin], step], 2), 'Snow melt:', - np.round(self.bin_meltsnow[glac_idx_t0[0], step], 2), + np.round(self.bin_meltsnow[glac_idx_t0[nbin], step], 2), 'Rain:', - np.round(self.bin_prec[glac_idx_t0[0], step], 2), + np.round(self.bin_prec[glac_idx_t0[nbin], step], 2), + 'Rfrz:', + np.round(self.bin_refreeze[gidx, step], 2), ) - # Refreeze [m w.e.] - # refreeze cannot exceed rain and melt (snow & glacier melt) - self.bin_refreeze[:, step] = self.bin_meltsnow[:, step] + self.bin_prec[:, step] - # refreeze cannot exceed snow depth - self.bin_refreeze[ - self.bin_refreeze[:, step] > self.bin_snowpack[:, step], - step, - ] = self.bin_snowpack[ - self.bin_refreeze[:, step] > self.bin_snowpack[:, step], - step, - ] - # refreeze cannot exceed refreeze potential - self.bin_refreeze[self.bin_refreeze[:, step] > refreeze_potential, step] = refreeze_potential[ - self.bin_refreeze[:, step] > refreeze_potential - ] - self.bin_refreeze[ - abs(self.bin_refreeze[:, step]) < pygem_prms['constants']['tolerance'], - step, - ] = 0 - # update refreeze potential - refreeze_potential -= self.bin_refreeze[:, step] - refreeze_potential[abs(refreeze_potential) < pygem_prms['constants']['tolerance']] = 0 - - # SNOWPACK REMAINING [m w.e.] - self.snowpack_remaining[:, step] = self.bin_snowpack[:, step] - self.bin_meltsnow[:, step] - self.snowpack_remaining[ - abs(self.snowpack_remaining[:, step]) < pygem_prms['constants']['tolerance'], + elif pygem_prms['mb']['option_refreezing'] == 'Woodward': + # Refreeze based on annual air temperature (Woodward etal. 1997) + # R(m) = (-0.69 * Tair + 0.0096) * 1 m / 100 cm + # calculate annually and place potential refreeze in user defined month + if step == t_start: + bin_temp_annual = annualweightedmean_array( + self.bin_temp[:, t_start : t_stop + 1], + self.dates_table.iloc[t_start : t_stop + 1, :], + ) + bin_refreezepotential_annual = (-0.69 * bin_temp_annual + 0.0096) / 100 + # Remove negative refreezing values + bin_refreezepotential_annual[bin_refreezepotential_annual < 0] = 0 + self.bin_refreezepotential[:, step] = bin_refreezepotential_annual + # Reset refreeze potential every year + if self.bin_refreezepotential[:, step].max() > 0: + refreeze_potential = self.bin_refreezepotential[:, step] + + if self.debug_refreeze: + print( + self.dates_table.loc[step, 'date'], + 'Refreeze potential:', + np.round(refreeze_potential[glac_idx_t0[50]], 3), + 'Snow depth:', + np.round(self.bin_snowpack[glac_idx_t0[50], step], 2), + 'Snow melt:', + np.round(self.bin_meltsnow[glac_idx_t0[50], step], 2), + 'Rain:', + np.round(self.bin_prec[glac_idx_t0[50], step], 2), + ) + + # Refreeze [m w.e.] + # refreeze cannot exceed rain and melt (snow & glacier melt) + self.bin_refreeze[:, step] = self.bin_meltsnow[:, step] + self.bin_prec[:, step] + # refreeze cannot exceed snow depth + self.bin_refreeze[self.bin_refreeze[:, step] > self.bin_snowpack[:, step], step] = ( + self.bin_snowpack[self.bin_refreeze[:, step] > self.bin_snowpack[:, step], step] + ) + # refreeze cannot exceed refreeze potential + self.bin_refreeze[self.bin_refreeze[:, step] > refreeze_potential, step] = refreeze_potential[ + self.bin_refreeze[:, step] > refreeze_potential + ] + self.bin_refreeze[ + abs(self.bin_refreeze[:, step]) < pygem_prms['constants']['tolerance'], step, ] = 0 + # update refreeze potential + refreeze_potential -= self.bin_refreeze[:, step] + refreeze_potential[abs(refreeze_potential) < pygem_prms['constants']['tolerance']] = 0 + + # SNOWPACK REMAINING [m w.e.] + self.snowpack_remaining[:, step] = self.bin_snowpack[:, step] - self.bin_meltsnow[:, step] + self.snowpack_remaining[ + abs(self.snowpack_remaining[:, step]) < pygem_prms['constants']['tolerance'], + step, + ] = 0 + + # Record values + self.glac_bin_melt[glac_idx_t0, step] = self.bin_melt[glac_idx_t0, step] + self.glac_bin_refreeze[glac_idx_t0, step] = self.bin_refreeze[glac_idx_t0, step] + self.glac_bin_snowpack[glac_idx_t0, step] = self.bin_snowpack[glac_idx_t0, step] + # CLIMATIC MASS BALANCE [m w.e.] + self.glac_bin_massbalclim[glac_idx_t0, step] = ( + self.bin_acc[glac_idx_t0, step] + + self.glac_bin_refreeze[glac_idx_t0, step] + - self.glac_bin_melt[glac_idx_t0, step] + ) - # Record values - self.glac_bin_melt[glac_idx_t0, step] = self.bin_melt[glac_idx_t0, step] - self.glac_bin_refreeze[glac_idx_t0, step] = self.bin_refreeze[glac_idx_t0, step] - self.glac_bin_snowpack[glac_idx_t0, step] = self.bin_snowpack[glac_idx_t0, step] - # CLIMATIC MASS BALANCE [m w.e.] - self.glac_bin_massbalclim[glac_idx_t0, step] = ( - self.bin_acc[glac_idx_t0, step] - + self.glac_bin_refreeze[glac_idx_t0, step] - - self.glac_bin_melt[glac_idx_t0, step] + # OFF-GLACIER ACCUMULATION, MELT, REFREEZE, AND SNOWPACK + if option_areaconstant == False: + # precipitation, refreeze, and snowpack are the same both on- and off-glacier + self.offglac_bin_prec[offglac_idx, step] = self.bin_prec[offglac_idx, step] + self.offglac_bin_refreeze[offglac_idx, step] = self.bin_refreeze[offglac_idx, step] + self.offglac_bin_snowpack[offglac_idx, step] = self.bin_snowpack[offglac_idx, step] + # Off-glacier melt includes both snow melt and melting of refreezing + # (this is not an issue on-glacier because energy remaining melts underlying snow/ice) + # melt of refreezing (assumed to be snow) + self.offglac_meltrefreeze = self.surfacetype_ddf_dict[2] * melt_energy_available + # melt of refreezing cannot exceed refreezing + self.offglac_meltrefreeze[self.offglac_meltrefreeze > self.bin_refreeze[:, step]] = ( + self.bin_refreeze[:, step][self.offglac_meltrefreeze > self.bin_refreeze[:, step]] + ) + # off-glacier melt = snow melt + refreezing melt + self.offglac_bin_melt[offglac_idx, step] = ( + self.bin_meltsnow[offglac_idx, step] + self.offglac_meltrefreeze[offglac_idx] ) - # OFF-GLACIER ACCUMULATION, MELT, REFREEZE, AND SNOWPACK - if option_areaconstant == False: - # precipitation, refreeze, and snowpack are the same both on- and off-glacier - self.offglac_bin_prec[offglac_idx, step] = self.bin_prec[offglac_idx, step] - self.offglac_bin_refreeze[offglac_idx, step] = self.bin_refreeze[offglac_idx, step] - self.offglac_bin_snowpack[offglac_idx, step] = self.bin_snowpack[offglac_idx, step] - # Off-glacier melt includes both snow melt and melting of refreezing - # (this is not an issue on-glacier because energy remaining melts underlying snow/ice) - # melt of refreezing (assumed to be snow) - self.offglac_meltrefreeze = self.surfacetype_ddf_dict[2] * melt_energy_available - # melt of refreezing cannot exceed refreezing - self.offglac_meltrefreeze[self.offglac_meltrefreeze > self.bin_refreeze[:, step]] = ( - self.bin_refreeze[:, step][self.offglac_meltrefreeze > self.bin_refreeze[:, step]] - ) - # off-glacier melt = snow melt + refreezing melt - self.offglac_bin_melt[offglac_idx, step] = ( - self.bin_meltsnow[offglac_idx, step] + self.offglac_meltrefreeze[offglac_idx] - ) - - # ===== RETURN TO ANNUAL LOOP ===== - # SURFACE TYPE (-) - # Annual climatic mass balance [m w.e.] used to determine the surface type - self.glac_bin_massbalclim_annual[:, year_idx] = self.glac_bin_massbalclim[ - :, year_start_month_idx:year_stop_month_idx - ].sum(1) - # Update surface type for each bin - self.surfacetype, firnline_idx = self._surfacetypebinsannual( - self.surfacetype, self.glac_bin_massbalclim_annual, year_idx - ) - # Record binned glacier area - self.glac_bin_area_annual[:, year_idx] = glacier_area_t0 - # Store glacier-wide results - self._convert_glacwide_results( - year_idx, - year_start_month_idx, - year_stop_month_idx, - glacier_area_t0, - heights, - fls=fls, - fl_id=fl_id, - option_areaconstant=option_areaconstant, - ) + # ===== RETURN TO ANNUAL LOOP ===== + # SURFACE TYPE (-) + # Annual climatic mass balance [m w.e.] used to determine the surface type + self.glac_bin_massbalclim_annual[:, year_idx] = self.glac_bin_massbalclim[:, t_start : t_stop + 1].sum(1) + # Update surface type for each bin + self.surfacetype, firnline_idx = self._surfacetypebinsannual( + self.surfacetype, self.glac_bin_massbalclim_annual, year_idx + ) + # Record binned glacier area + self.glac_bin_area_annual[:, year_idx] = glacier_area_t0 + # Store glacier-wide results + self._convert_glacwide_results( + year_idx, + t_start, + t_stop, + glacier_area_t0, + heights, + fls=fls, + fl_id=fl_id, + option_areaconstant=option_areaconstant, + ) - # Mass balance for each bin [m ice per second] - seconds_in_year = self.dayspermonth[year_start_month_idx:year_stop_month_idx].sum() * 24 * 3600 + # Calculate mass balance rate for each bin (convert from [m w.e. yr^-1] to [m ice s^-1]) + seconds_in_year = self.days_in_step[t_start : t_stop + 1].sum() * 24 * 3600 mb = ( - self.glac_bin_massbalclim[:, year_start_month_idx:year_stop_month_idx].sum(1) + self.glac_bin_massbalclim[:, t_start : t_stop + 1].sum(1) * pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice'] / seconds_in_year @@ -799,7 +796,6 @@ def get_annual_mb( mb_filled[(mb_filled == 0) & (heights < height_max)] = mb_min + mb_grad * ( height_min - heights[(mb_filled == 0) & (heights < height_max)] ) - elif len(glac_idx_t0) >= 1 and len(glac_idx_t0) <= 3 and mb.max() <= 0: mb_min = np.min(mb[glac_idx_t0]) height_max = np.max(heights[glac_idx_t0]) @@ -811,8 +807,8 @@ def get_annual_mb( def _convert_glacwide_results( self, year_idx, - year_start_month_idx, - year_stop_month_idx, + t_start, + t_stop, glacier_area, heights, fls=None, @@ -838,7 +834,8 @@ def _convert_glacwide_results( """ # Glacier area glac_idx = glacier_area.nonzero()[0] - glacier_area_monthly = glacier_area[:, np.newaxis].repeat(12, axis=1) + # repeat glacier area for all time steps in year + glacier_area_steps = glacier_area[:, np.newaxis].repeat((t_stop + 1) - t_start, axis=1) # Check if need to adjust for complete removal of the glacier # - needed for accurate runoff calcs and accurate mass balance components @@ -855,7 +852,7 @@ def _convert_glacwide_results( ) # Check annual climatic mass balance (mwea) mb_mwea = ( - glacier_area * self.glac_bin_massbalclim[:, year_start_month_idx:year_stop_month_idx].sum(1) + glacier_area * self.glac_bin_massbalclim[:, t_start : t_stop + 1].sum(1) ).sum() / glacier_area.sum() else: mb_max_loss = 0 @@ -876,82 +873,77 @@ def _convert_glacwide_results( # Glacier-wide area (m2) self.glac_wide_area_annual[year_idx] = glacier_area.sum() # Glacier-wide temperature (degC) - self.glac_wide_temp[year_start_month_idx:year_stop_month_idx] = ( - self.bin_temp[:, year_start_month_idx:year_stop_month_idx][glac_idx] * glacier_area_monthly[glac_idx] + self.glac_wide_temp[t_start : t_stop + 1] = ( + self.bin_temp[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) / glacier_area.sum() # Glacier-wide precipitation (m3) - self.glac_wide_prec[year_start_month_idx:year_stop_month_idx] = ( - self.bin_prec[:, year_start_month_idx:year_stop_month_idx][glac_idx] * glacier_area_monthly[glac_idx] + self.glac_wide_prec[t_start : t_stop + 1] = ( + self.bin_prec[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) # Glacier-wide accumulation (m3 w.e.) - self.glac_wide_acc[year_start_month_idx:year_stop_month_idx] = ( - self.bin_acc[:, year_start_month_idx:year_stop_month_idx][glac_idx] * glacier_area_monthly[glac_idx] + self.glac_wide_acc[t_start : t_stop + 1] = ( + self.bin_acc[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) # Glacier-wide refreeze (m3 w.e.) - self.glac_wide_refreeze[year_start_month_idx:year_stop_month_idx] = ( - self.glac_bin_refreeze[:, year_start_month_idx:year_stop_month_idx][glac_idx] - * glacier_area_monthly[glac_idx] + self.glac_wide_refreeze[t_start : t_stop + 1] = ( + self.glac_bin_refreeze[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) # Glacier-wide melt (m3 w.e.) - self.glac_wide_melt[year_start_month_idx:year_stop_month_idx] = ( - self.glac_bin_melt[:, year_start_month_idx:year_stop_month_idx][glac_idx] - * glacier_area_monthly[glac_idx] + self.glac_wide_melt[t_start : t_stop + 1] = ( + self.glac_bin_melt[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) # Glacier-wide total mass balance (m3 w.e.) - self.glac_wide_massbaltotal[year_start_month_idx:year_stop_month_idx] = ( - self.glac_wide_acc[year_start_month_idx:year_stop_month_idx] - + self.glac_wide_refreeze[year_start_month_idx:year_stop_month_idx] - - self.glac_wide_melt[year_start_month_idx:year_stop_month_idx] - - self.glac_wide_frontalablation[year_start_month_idx:year_stop_month_idx] + self.glac_wide_massbaltotal[t_start : t_stop + 1] = ( + self.glac_wide_acc[t_start : t_stop + 1] + + self.glac_wide_refreeze[t_start : t_stop + 1] + - self.glac_wide_melt[t_start : t_stop + 1] + - self.glac_wide_frontalablation[t_start : t_stop + 1] ) # If mass loss more negative than glacier mass, reduce melt so glacier completely melts (no excess) if icethickness_t0 is not None and mb_mwea < mb_max_loss: - melt_yr_raw = self.glac_wide_melt[year_start_month_idx:year_stop_month_idx].sum() + melt_yr_raw = self.glac_wide_melt[t_start : t_stop + 1].sum() melt_yr_max = ( self.glac_wide_volume_annual[year_idx] * pygem_prms['constants']['density_ice'] / pygem_prms['constants']['density_water'] - + self.glac_wide_acc[year_start_month_idx:year_stop_month_idx].sum() - + self.glac_wide_refreeze[year_start_month_idx:year_stop_month_idx].sum() + + self.glac_wide_acc[t_start : t_stop + 1].sum() + + self.glac_wide_refreeze[t_start : t_stop + 1].sum() ) melt_frac = melt_yr_max / melt_yr_raw # Update glacier-wide melt (m3 w.e.) - self.glac_wide_melt[year_start_month_idx:year_stop_month_idx] = ( - self.glac_wide_melt[year_start_month_idx:year_stop_month_idx] * melt_frac - ) + self.glac_wide_melt[t_start : t_stop + 1] = self.glac_wide_melt[t_start : t_stop + 1] * melt_frac # Glacier-wide runoff (m3) - self.glac_wide_runoff[year_start_month_idx:year_stop_month_idx] = ( - self.glac_wide_prec[year_start_month_idx:year_stop_month_idx] - + self.glac_wide_melt[year_start_month_idx:year_stop_month_idx] - - self.glac_wide_refreeze[year_start_month_idx:year_stop_month_idx] + self.glac_wide_runoff[t_start : t_stop + 1] = ( + self.glac_wide_prec[t_start : t_stop + 1] + + self.glac_wide_melt[t_start : t_stop + 1] + - self.glac_wide_refreeze[t_start : t_stop + 1] ) + # Snow line altitude (m a.s.l.) - heights_monthly = heights[:, np.newaxis].repeat(12, axis=1) - snow_mask = np.zeros(heights_monthly.shape) - snow_mask[self.glac_bin_snowpack[:, year_start_month_idx:year_stop_month_idx] > 0] = 1 - heights_monthly_wsnow = heights_monthly * snow_mask - heights_monthly_wsnow[heights_monthly_wsnow == 0] = np.nan + heights_steps = heights[:, np.newaxis].repeat((t_stop + 1) - t_start, axis=1) + snow_mask = np.zeros(heights_steps.shape) + snow_mask[self.glac_bin_snowpack[:, t_start : t_stop + 1] > 0] = 1 + heights_steps_wsnow = heights_steps * snow_mask + heights_steps_wsnow[heights_steps_wsnow == 0] = np.nan heights_change = np.zeros(heights.shape) heights_change[0:-1] = heights[0:-1] - heights[1:] try: - snowline_idx = np.nanargmin(heights_monthly_wsnow, axis=0) - self.glac_wide_snowline[year_start_month_idx:year_stop_month_idx] = ( - heights[snowline_idx] - heights_change[snowline_idx] / 2 - ) + snowline_idx = np.nanargmin(heights_steps_wsnow, axis=0) + self.glac_wide_snowline[t_start : t_stop + 1] = heights[snowline_idx] - heights_change[snowline_idx] / 2 except: - snowline_idx = np.zeros((heights_monthly_wsnow.shape[1])).astype(int) + snowline_idx = np.zeros((heights_steps_wsnow.shape[1])).astype(int) snowline_idx_nan = [] - for ncol in range(heights_monthly_wsnow.shape[1]): - if ~np.isnan(heights_monthly_wsnow[:, ncol]).all(): - snowline_idx[ncol] = np.nanargmin(heights_monthly_wsnow[:, ncol]) + for ncol in range(heights_steps_wsnow.shape[1]): + if ~np.isnan(heights_steps_wsnow[:, ncol]).all(): + snowline_idx[ncol] = np.nanargmin(heights_steps_wsnow[:, ncol]) else: snowline_idx_nan.append(ncol) heights_manual = heights[snowline_idx] - heights_change[snowline_idx] / 2 heights_manual[snowline_idx_nan] = np.nan # this line below causes a potential All-NaN slice encountered issue at some time steps - self.glac_wide_snowline[year_start_month_idx:year_stop_month_idx] = heights_manual + self.glac_wide_snowline[t_start : t_stop + 1] = heights_manual # Equilibrium line altitude (m a.s.l.) ela_mask = np.zeros(heights.shape) @@ -967,33 +959,31 @@ def _convert_glacwide_results( # ===== Off-glacier ==== offglac_idx = np.where(self.offglac_bin_area_annual[:, year_idx] > 0)[0] if option_areaconstant == False and len(offglac_idx) > 0: - offglacier_area_monthly = self.offglac_bin_area_annual[:, year_idx][:, np.newaxis].repeat(12, axis=1) + offglacier_area_steps = self.offglac_bin_area_annual[:, year_idx][:, np.newaxis].repeat( + (t_stop + 1) - t_start, axis=1 + ) # Off-glacier precipitation (m3) - self.offglac_wide_prec[year_start_month_idx:year_stop_month_idx] = ( - self.bin_prec[:, year_start_month_idx:year_stop_month_idx][offglac_idx] - * offglacier_area_monthly[offglac_idx] + self.offglac_wide_prec[t_start : t_stop + 1] = ( + self.bin_prec[:, t_start : t_stop + 1][offglac_idx] * offglacier_area_steps[offglac_idx] ).sum(0) # Off-glacier melt (m3 w.e.) - self.offglac_wide_melt[year_start_month_idx:year_stop_month_idx] = ( - self.offglac_bin_melt[:, year_start_month_idx:year_stop_month_idx][offglac_idx] - * offglacier_area_monthly[offglac_idx] + self.offglac_wide_melt[t_start : t_stop + 1] = ( + self.offglac_bin_melt[:, t_start : t_stop + 1][offglac_idx] * offglacier_area_steps[offglac_idx] ).sum(0) # Off-glacier refreeze (m3 w.e.) - self.offglac_wide_refreeze[year_start_month_idx:year_stop_month_idx] = ( - self.offglac_bin_refreeze[:, year_start_month_idx:year_stop_month_idx][offglac_idx] - * offglacier_area_monthly[offglac_idx] + self.offglac_wide_refreeze[t_start : t_stop + 1] = ( + self.offglac_bin_refreeze[:, t_start : t_stop + 1][offglac_idx] * offglacier_area_steps[offglac_idx] ).sum(0) # Off-glacier runoff (m3) - self.offglac_wide_runoff[year_start_month_idx:year_stop_month_idx] = ( - self.offglac_wide_prec[year_start_month_idx:year_stop_month_idx] - + self.offglac_wide_melt[year_start_month_idx:year_stop_month_idx] - - self.offglac_wide_refreeze[year_start_month_idx:year_stop_month_idx] + self.offglac_wide_runoff[t_start : t_stop + 1] = ( + self.offglac_wide_prec[t_start : t_stop + 1] + + self.offglac_wide_melt[t_start : t_stop + 1] + - self.offglac_wide_refreeze[t_start : t_stop + 1] ) # Off-glacier snowpack (m3 w.e.) - self.offglac_wide_snowpack[year_start_month_idx:year_stop_month_idx] = ( - self.offglac_bin_snowpack[:, year_start_month_idx:year_stop_month_idx][offglac_idx] - * offglacier_area_monthly[offglac_idx] + self.offglac_wide_snowpack[t_start : t_stop + 1] = ( + self.offglac_bin_snowpack[:, t_start : t_stop + 1][offglac_idx] * offglacier_area_steps[offglac_idx] ).sum(0) def ensure_mass_conservation(self, diag): @@ -1010,12 +1000,28 @@ def ensure_mass_conservation(self, diag): total volume change and therefore do not impose limitations like this because they do not estimate the flux divergence. As a result, they may systematically overestimate mass loss compared to OGGM's dynamical model. """ + years = list(np.unique(self.dates_table['year'])) + + # Compute annual volume change and melt values needed for adjustments + vol_change_annual_mbmod = np.zeros(len(years)) + vol_change_annual_mbmod_melt = np.zeros(len(years)) + for nyear, year in enumerate(years): + # get time step indices - note indexing should be [t_start:t_stop+1] to include final step in year + t_start, t_stop = self.get_step_inds(year) + + vol_change_annual_mbmod[nyear] = ( + self.glac_wide_massbaltotal[t_start : t_stop + 1].sum() + * pygem_prms['constants']['density_water'] + / pygem_prms['constants']['density_ice'] + ) + + vol_change_annual_mbmod_melt[nyear] = ( + self.glac_wide_melt[t_start : t_stop + 1].sum() + * pygem_prms['constants']['density_water'] + / pygem_prms['constants']['density_ice'] + ) + # Compute difference between volume change - vol_change_annual_mbmod = ( - self.glac_wide_massbaltotal.reshape(-1, 12).sum(1) - * pygem_prms['constants']['density_water'] - / pygem_prms['constants']['density_ice'] - ) vol_change_annual_diag = np.zeros(vol_change_annual_mbmod.shape) vol_change_annual_diag[0 : diag.volume_m3.values[1:].shape[0]] = ( diag.volume_m3.values[1:] - diag.volume_m3.values[:-1] @@ -1023,11 +1029,6 @@ def ensure_mass_conservation(self, diag): vol_change_annual_dif = vol_change_annual_diag - vol_change_annual_mbmod # Reduce glacier melt by the difference - vol_change_annual_mbmod_melt = ( - self.glac_wide_melt.reshape(-1, 12).sum(1) - * pygem_prms['constants']['density_water'] - / pygem_prms['constants']['density_ice'] - ) vol_change_annual_melt_reduction = np.zeros(vol_change_annual_mbmod.shape) chg_idx = vol_change_annual_mbmod.nonzero()[0] chg_idx_posmbmod = vol_change_annual_mbmod_melt.nonzero()[0] @@ -1037,7 +1038,11 @@ def ensure_mass_conservation(self, diag): 1 - vol_change_annual_dif[chg_idx_melt] / vol_change_annual_mbmod_melt[chg_idx_melt] ) - vol_change_annual_melt_reduction_monthly = np.repeat(vol_change_annual_melt_reduction, 12) + vol_change_annual_melt_reduction_monthly = np.zeros(self.dates_table.shape[0]) + for nyear, year in enumerate(years): + # get time step indices - note indexing should be [t_start:t_stop+1] to include final step in year + t_start, t_stop = self.get_step_inds(year) + vol_change_annual_melt_reduction_monthly[t_start : t_stop + 1] = vol_change_annual_melt_reduction[nyear] # Glacier-wide melt (m3 w.e.) self.glac_wide_melt = self.glac_wide_melt * vol_change_annual_melt_reduction_monthly diff --git a/pygem/output.py b/pygem/output.py index c76005f2..f1cfa7f8 100644 --- a/pygem/output.py +++ b/pygem/output.py @@ -53,6 +53,8 @@ class single_glacier: DataFrame containing metadata and characteristics of the glacier from the Randolph Glacier Inventory. dates_table : pd.DataFrame DataFrame containing the time series of dates associated with the model output. + timestep : str + The time step resolution ('monthly' or 'daily') sim_climate_name : str Name of the General Circulation Model (GCM) used for climate forcing. sim_climate_scenario : str @@ -79,6 +81,7 @@ class single_glacier: glacier_rgi_table: pd.DataFrame dates_table: pd.DataFrame + timestep: str sim_climate_name: str sim_climate_scenario: str realization: str @@ -90,6 +93,7 @@ class single_glacier: sim_endyear: int option_calibration: str option_bias_adjustment: str + option_dynamics: str extra_vars: bool = False def __post_init__(self): @@ -176,8 +180,7 @@ def _set_time_vals(self): self.annual_columns = np.unique(self.dates_table['year'].values)[0 : int(self.dates_table.shape[0] / 12)] elif pygem_prms['climate']['sim_wateryear'] == 'custom': self.year_type = 'custom year' - self.time_values = self.dates_table['date'].values.tolist() - self.time_values = [cftime.DatetimeNoLeap(x.year, x.month, x.day) for x in self.time_values] + self.time_values = [cftime.DatetimeGregorian(x.year, x.month, x.day) for x in self.dates_table['date']] # append additional year to self.year_values to account for mass and area at end of period self.year_values = self.annual_columns self.year_values = np.concatenate((self.year_values, np.array([self.annual_columns[-1] + 1]))) @@ -196,6 +199,8 @@ def _model_params_record(self): self.mdl_params_dict['sim_climate_scenario'] = self.sim_climate_scenario self.mdl_params_dict['option_calibration'] = self.option_calibration self.mdl_params_dict['option_bias_adjustment'] = self.option_bias_adjustment + self.mdl_params_dict['option_dynamics'] = self.option_dynamics + self.mdl_params_dict['timestep'] = self.timestep # record manually defined modelprms if calibration option is None if not self.option_calibration: self._update_modelparams_record() @@ -349,13 +354,13 @@ def _set_outdir(self): def _update_dicts(self): """Update coordinate and attribute dictionaries specific to glacierwide_stats outputs""" - self.output_coords_dict['glac_runoff_monthly'] = collections.OrderedDict( + self.output_coords_dict['glac_runoff'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_runoff_monthly'] = { + self.output_attrs_dict['glac_runoff'] = { 'long_name': 'glacier-wide runoff', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'runoff from the glacier terminus, which moves over time', } self.output_coords_dict['glac_area_annual'] = collections.OrderedDict( @@ -394,25 +399,25 @@ def _update_dicts(self): 'temporal_resolution': 'annual', 'comment': 'equilibrium line altitude is the elevation where the climatic mass balance is zero', } - self.output_coords_dict['offglac_runoff_monthly'] = collections.OrderedDict( + self.output_coords_dict['offglac_runoff'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_runoff_monthly'] = { + self.output_attrs_dict['offglac_runoff'] = { 'long_name': 'off-glacier-wide runoff', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'off-glacier runoff from area where glacier no longer exists', } # if nsims > 1, store median-absolute deviation metrics if self.nsims > 1: - self.output_coords_dict['glac_runoff_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['glac_runoff_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_runoff_monthly_mad'] = { + self.output_attrs_dict['glac_runoff_mad'] = { 'long_name': 'glacier-wide runoff median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'runoff from the glacier terminus, which moves over time', } self.output_coords_dict['glac_area_annual_mad'] = collections.OrderedDict( @@ -451,93 +456,93 @@ def _update_dicts(self): 'temporal_resolution': 'annual', 'comment': 'equilibrium line altitude is the elevation where the climatic mass balance is zero', } - self.output_coords_dict['offglac_runoff_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['offglac_runoff_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_runoff_monthly_mad'] = { + self.output_attrs_dict['offglac_runoff_mad'] = { 'long_name': 'off-glacier-wide runoff median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'off-glacier runoff from area where glacier no longer exists', } # optionally store extra variables if self.extra_vars: - self.output_coords_dict['glac_prec_monthly'] = collections.OrderedDict( + self.output_coords_dict['glac_prec'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_prec_monthly'] = { + self.output_attrs_dict['glac_prec'] = { 'long_name': 'glacier-wide precipitation (liquid)', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'only the liquid precipitation, solid precipitation excluded', } - self.output_coords_dict['glac_temp_monthly'] = collections.OrderedDict( + self.output_coords_dict['glac_temp'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_temp_monthly'] = { + self.output_attrs_dict['glac_temp'] = { 'standard_name': 'air_temperature', 'long_name': 'glacier-wide mean air temperature', 'units': 'K', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': ( 'each elevation bin is weighted equally to compute the mean temperature, and ' 'bins where the glacier no longer exists due to retreat have been removed' ), } - self.output_coords_dict['glac_acc_monthly'] = collections.OrderedDict( + self.output_coords_dict['glac_acc'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_acc_monthly'] = { + self.output_attrs_dict['glac_acc'] = { 'long_name': 'glacier-wide accumulation, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'only the solid precipitation', } - self.output_coords_dict['glac_refreeze_monthly'] = collections.OrderedDict( + self.output_coords_dict['glac_refreeze'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_refreeze_monthly'] = { + self.output_attrs_dict['glac_refreeze'] = { 'long_name': 'glacier-wide refreeze, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, } - self.output_coords_dict['glac_melt_monthly'] = collections.OrderedDict( + self.output_coords_dict['glac_melt'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_melt_monthly'] = { + self.output_attrs_dict['glac_melt'] = { 'long_name': 'glacier-wide melt, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, } - self.output_coords_dict['glac_frontalablation_monthly'] = collections.OrderedDict( + self.output_coords_dict['glac_frontalablation'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_frontalablation_monthly'] = { + self.output_attrs_dict['glac_frontalablation'] = { 'long_name': 'glacier-wide frontal ablation, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': ( 'mass losses from calving, subaerial frontal melting, sublimation above the ' 'waterline and subaqueous frontal melting below the waterline; positive values indicate mass lost like melt' ), } - self.output_coords_dict['glac_massbaltotal_monthly'] = collections.OrderedDict( + self.output_coords_dict['glac_massbaltotal'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_massbaltotal_monthly'] = { + self.output_attrs_dict['glac_massbaltotal'] = { 'long_name': 'glacier-wide total mass balance, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'total mass balance is the sum of the climatic mass balance and frontal ablation', } - self.output_coords_dict['glac_snowline_monthly'] = collections.OrderedDict( + self.output_coords_dict['glac_snowline'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_snowline_monthly'] = { + self.output_attrs_dict['glac_snowline'] = { 'long_name': 'transient snowline altitude above mean sea level', 'units': 'm', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'transient snowline is altitude separating snow from ice/firn', } self.output_coords_dict['glac_mass_change_ignored_annual'] = collections.OrderedDict( @@ -549,119 +554,119 @@ def _update_dicts(self): 'temporal_resolution': 'annual', 'comment': 'glacier mass change ignored due to flux divergence', } - self.output_coords_dict['offglac_prec_monthly'] = collections.OrderedDict( + self.output_coords_dict['offglac_prec'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_prec_monthly'] = { + self.output_attrs_dict['offglac_prec'] = { 'long_name': 'off-glacier-wide precipitation (liquid)', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'only the liquid precipitation, solid precipitation excluded', } - self.output_coords_dict['offglac_refreeze_monthly'] = collections.OrderedDict( + self.output_coords_dict['offglac_refreeze'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_refreeze_monthly'] = { + self.output_attrs_dict['offglac_refreeze'] = { 'long_name': 'off-glacier-wide refreeze, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, } - self.output_coords_dict['offglac_melt_monthly'] = collections.OrderedDict( + self.output_coords_dict['offglac_melt'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_melt_monthly'] = { + self.output_attrs_dict['offglac_melt'] = { 'long_name': 'off-glacier-wide melt, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'only melt of snow and refreeze since off-glacier', } - self.output_coords_dict['offglac_snowpack_monthly'] = collections.OrderedDict( + self.output_coords_dict['offglac_snowpack'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_snowpack_monthly'] = { + self.output_attrs_dict['offglac_snowpack'] = { 'long_name': 'off-glacier-wide snowpack, in water equivalent', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'snow remaining accounting for new accumulation, melt, and refreeze', } # if nsims > 1, store median-absolute deviation metrics if self.nsims > 1: - self.output_coords_dict['glac_prec_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['glac_prec_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_prec_monthly_mad'] = { + self.output_attrs_dict['glac_prec_mad'] = { 'long_name': 'glacier-wide precipitation (liquid) median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'only the liquid precipitation, solid precipitation excluded', } - self.output_coords_dict['glac_temp_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['glac_temp_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_temp_monthly_mad'] = { + self.output_attrs_dict['glac_temp_mad'] = { 'standard_name': 'air_temperature', 'long_name': 'glacier-wide mean air temperature median absolute deviation', 'units': 'K', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': ( 'each elevation bin is weighted equally to compute the mean temperature, and ' 'bins where the glacier no longer exists due to retreat have been removed' ), } - self.output_coords_dict['glac_acc_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['glac_acc_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_acc_monthly_mad'] = { + self.output_attrs_dict['glac_acc_mad'] = { 'long_name': 'glacier-wide accumulation, in water equivalent, median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'only the solid precipitation', } - self.output_coords_dict['glac_refreeze_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['glac_refreeze_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_refreeze_monthly_mad'] = { + self.output_attrs_dict['glac_refreeze_mad'] = { 'long_name': 'glacier-wide refreeze, in water equivalent, median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, } - self.output_coords_dict['glac_melt_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['glac_melt_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_melt_monthly_mad'] = { + self.output_attrs_dict['glac_melt_mad'] = { 'long_name': 'glacier-wide melt, in water equivalent, median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, } - self.output_coords_dict['glac_frontalablation_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['glac_frontalablation_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_frontalablation_monthly_mad'] = { + self.output_attrs_dict['glac_frontalablation_mad'] = { 'long_name': 'glacier-wide frontal ablation, in water equivalent, median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': ( 'mass losses from calving, subaerial frontal melting, sublimation above the ' 'waterline and subaqueous frontal melting below the waterline' ), } - self.output_coords_dict['glac_massbaltotal_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['glac_massbaltotal_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_massbaltotal_monthly_mad'] = { + self.output_attrs_dict['glac_massbaltotal_mad'] = { 'long_name': 'glacier-wide total mass balance, in water equivalent, median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'total mass balance is the sum of the climatic mass balance and frontal ablation', } - self.output_coords_dict['glac_snowline_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['glac_snowline_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['glac_snowline_monthly_mad'] = { + self.output_attrs_dict['glac_snowline_mad'] = { 'long_name': 'transient snowline above mean sea level median absolute deviation', 'units': 'm', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'transient snowline is altitude separating snow from ice/firn', } self.output_coords_dict['glac_mass_change_ignored_annual_mad'] = collections.OrderedDict( @@ -673,39 +678,39 @@ def _update_dicts(self): 'temporal_resolution': 'annual', 'comment': 'glacier mass change ignored due to flux divergence', } - self.output_coords_dict['offglac_prec_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['offglac_prec_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_prec_monthly_mad'] = { + self.output_attrs_dict['offglac_prec_mad'] = { 'long_name': 'off-glacier-wide precipitation (liquid) median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'only the liquid precipitation, solid precipitation excluded', } - self.output_coords_dict['offglac_refreeze_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['offglac_refreeze_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_refreeze_monthly_mad'] = { + self.output_attrs_dict['offglac_refreeze_mad'] = { 'long_name': 'off-glacier-wide refreeze, in water equivalent, median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, } - self.output_coords_dict['offglac_melt_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['offglac_melt_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_melt_monthly_mad'] = { + self.output_attrs_dict['offglac_melt_mad'] = { 'long_name': 'off-glacier-wide melt, in water equivalent, median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'only melt of snow and refreeze since off-glacier', } - self.output_coords_dict['offglac_snowpack_monthly_mad'] = collections.OrderedDict( + self.output_coords_dict['offglac_snowpack_mad'] = collections.OrderedDict( [('glac', self.glac_values), ('time', self.time_values)] ) - self.output_attrs_dict['offglac_snowpack_monthly_mad'] = { + self.output_attrs_dict['offglac_snowpack_mad'] = { 'long_name': 'off-glacier-wide snowpack, in water equivalent, median absolute deviation', 'units': 'm3', - 'temporal_resolution': 'monthly', + 'temporal_resolution': self.timestep, 'comment': 'snow remaining accounting for new accumulation, melt, and refreeze', } @@ -824,60 +829,60 @@ def _update_dicts(self): 'comment': 'climatic mass balance is computed before dynamics so can theoretically exceed ice thickness', }, ) - self.output_coords_dict['bin_massbalclim_monthly'] = collections.OrderedDict( + self.output_coords_dict['bin_massbalclim'] = collections.OrderedDict( [ ('glac', self.glac_values), ('bin', self.bin_values), ('time', self.time_values), ] ) - self.output_attrs_dict['bin_massbalclim_monthly'] = { - 'long_name': 'binned monthly climatic mass balance, in water equivalent', + self.output_attrs_dict['bin_massbalclim'] = { + 'long_name': 'binned climatic mass balance, in water equivalent', 'units': 'm', - 'temporal_resolution': 'monthly', - 'comment': 'monthly climatic mass balance from the PyGEM mass balance module', + 'temporal_resolution': self.timestep, + 'comment': 'climatic mass balance from the PyGEM mass balance module', } # optionally store binned mass balance components if self.binned_components: - self.output_coords_dict['bin_accumulation_monthly'] = collections.OrderedDict( + self.output_coords_dict['bin_accumulation'] = collections.OrderedDict( [ ('glac', self.glac_values), ('bin', self.bin_values), ('time', self.time_values), ] ) - self.output_attrs_dict['bin_accumulation_monthly'] = { - 'long_name': 'binned monthly accumulation, in water equivalent', + self.output_attrs_dict['bin_accumulation'] = { + 'long_name': 'binned accumulation, in water equivalent', 'units': 'm', - 'temporal_resolution': 'monthly', - 'comment': 'monthly accumulation from the PyGEM mass balance module', + 'temporal_resolution': self.timestep, + 'comment': 'accumulation from the PyGEM mass balance module', } - self.output_coords_dict['bin_melt_monthly'] = collections.OrderedDict( + self.output_coords_dict['bin_melt'] = collections.OrderedDict( [ ('glac', self.glac_values), ('bin', self.bin_values), ('time', self.time_values), ] ) - self.output_attrs_dict['bin_melt_monthly'] = { - 'long_name': 'binned monthly melt, in water equivalent', + self.output_attrs_dict['bin_melt'] = { + 'long_name': 'binned melt, in water equivalent', 'units': 'm', - 'temporal_resolution': 'monthly', - 'comment': 'monthly melt from the PyGEM mass balance module', + 'temporal_resolution': self.timestep, + 'comment': 'melt from the PyGEM mass balance module', } - self.output_coords_dict['bin_refreeze_monthly'] = collections.OrderedDict( + self.output_coords_dict['bin_refreeze'] = collections.OrderedDict( [ ('glac', self.glac_values), ('bin', self.bin_values), ('time', self.time_values), ] ) - self.output_attrs_dict['bin_refreeze_monthly'] = { - 'long_name': 'binned monthly refreeze, in water equivalent', + self.output_attrs_dict['bin_refreeze'] = { + 'long_name': 'binned refreeze, in water equivalent', 'units': 'm', - 'temporal_resolution': 'monthly', - 'comment': 'monthly refreeze from the PyGEM mass balance module', + 'temporal_resolution': self.timestep, + 'comment': 'refreeze from the PyGEM mass balance module', } # if nsims > 1, store median-absolute deviation metrics diff --git a/pygem/plot/graphics.py b/pygem/plot/graphics.py index 74f0bbba..103725b3 100644 --- a/pygem/plot/graphics.py +++ b/pygem/plot/graphics.py @@ -18,7 +18,15 @@ from pygem.utils.stats import effective_n -def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): +def plot_modeloutput_section( + model=None, + ax=None, + title='', + lnlabel=None, + legendon=True, + lgdkwargs={'loc': 'upper right', 'fancybox': False, 'borderaxespad': 0, 'handlelength': 1}, + **kwargs, +): """Plots the result of the model output along the flowline. A paired down version of OGGMs graphics.plot_modeloutput_section() @@ -40,15 +48,12 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): ax = fig.add_axes([0.07, 0.08, 0.7, 0.84]) else: fig = plt.gcf() + # get n lines plotted on figure + nlines = len(plt.gca().get_lines()) - # Compute area histo - area = np.array([]) height = np.array([]) bed = np.array([]) for cls in fls: - a = cls.widths_m * cls.dx_meter * 1e-6 - a = np.where(cls.thick > 0, a, 0) - area = np.concatenate((area, a)) height = np.concatenate((height, cls.surface_h)) bed = np.concatenate((bed, cls.bed_h)) ylim = [bed.min(), height.max()] @@ -57,8 +62,11 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): cls = fls[-1] x = np.arange(cls.nx) * cls.dx * cls.map_dx - # Plot the bed - ax.plot(x, cls.bed_h, color='k', linewidth=2.5, label='Bed (Parab.)') + if nlines == 0: + if getattr(model, 'do_calving', False): + ax.hlines(model.water_level, x[0], x[-1], linestyles=':', label='Water level', color='C0') + # Plot the bed + ax.plot(x, cls.bed_h, color='k', linewidth=2.5, label='Bed (Parab.)') # Plot glacier t1 = cls.thick[:-2] @@ -77,7 +85,7 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): else: srfls = '-' - ax.plot(x, cls.surface_h, color=srfcolor, linewidth=2, ls=srfls, label='Glacier') + ax.plot(x, cls.surface_h, color=srfcolor, linewidth=2, ls=srfls, label=lnlabel) # Plot tributaries for i, inflow in zip(cls.inflow_indices, cls.inflows): @@ -99,16 +107,14 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): markeredgecolor='k', label='Tributary (inactive)', ) - if getattr(model, 'do_calving', False): - ax.hlines(model.water_level, x[0], x[-1], linestyles=':', color='C0') ax.set_ylim(ylim) - ax.spines['top'].set_color('none') ax.xaxis.set_ticks_position('bottom') ax.set_xlabel('Distance along flowline (m)') ax.set_ylabel('Altitude (m)') - + if legendon: + ax.legend(**lgdkwargs) # Title ax.set_title(title, loc='left') diff --git a/pygem/pygem_modelsetup.py b/pygem/pygem_modelsetup.py index 9524b878..18a7bb33 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -80,30 +80,30 @@ def datesmodelrun( # Select attributes of DateTimeIndex (dt.year, dt.month, and dt.daysinmonth) dates_table['year'] = dates_table['date'].dt.year dates_table['month'] = dates_table['date'].dt.month - dates_table['daysinmonth'] = dates_table['date'].dt.daysinmonth + dates_table['days_in_step'] = dates_table['date'].dt.daysinmonth dates_table['timestep'] = np.arange(len(dates_table['date'])) # Set date as index dates_table.set_index('timestep', inplace=True) # Remove leap year days if user selected this with option_leapyear if pygem_prms['time']['option_leapyear'] == 0: - mask1 = dates_table['daysinmonth'] == 29 - dates_table.loc[mask1, 'daysinmonth'] = 28 + mask1 = dates_table['days_in_step'] == 29 + dates_table.loc[mask1, 'days_in_step'] = 28 elif pygem_prms['time']['timestep'] == 'daily': # Automatically generate daily (freq = 'D') dates - dates_table = pd.DataFrame({'date': pd.date_range(startdate, enddate, freq='D')}) + dates_table = pd.DataFrame({'date': pd.date_range(startdate, enddate, freq='D', unit='s')}) # Extract attributes for dates_table dates_table['year'] = dates_table['date'].dt.year dates_table['month'] = dates_table['date'].dt.month dates_table['day'] = dates_table['date'].dt.day - dates_table['daysinmonth'] = dates_table['date'].dt.daysinmonth + dates_table['days_in_step'] = 1 dates_table['timestep'] = np.arange(len(dates_table['date'])) # Set date as index dates_table.set_index('timestep', inplace=True) # Remove leap year days if user selected this with option_leapyear if pygem_prms['time']['option_leapyear'] == 0: - # First, change 'daysinmonth' number - mask1 = dates_table['daysinmonth'] == 29 - dates_table.loc[mask1, 'daysinmonth'] = 28 + # First, change 'days_in_step' number + mask1 = dates_table['days_in_step'] == 29 + dates_table.loc[mask1, 'days_in_step'] = 28 # Next, remove the 29th days from the dates mask2 = (dates_table['month'] == 2) & (dates_table['day'] == 29) dates_table.drop(dates_table[mask2].index, inplace=True) @@ -135,52 +135,6 @@ def datesmodelrun( return dates_table -def daysinmonth(year, month): - """ - Return days in month based on the month and year - - Parameters - ---------- - year : str - month : str - - Returns - ------- - integer of the days in the month - """ - if year % 4 == 0: - daysinmonth_dict = { - 1: 31, - 2: 29, - 3: 31, - 4: 30, - 5: 31, - 6: 30, - 7: 31, - 8: 31, - 9: 30, - 10: 31, - 11: 30, - 12: 31, - } - else: - daysinmonth_dict = { - 1: 31, - 2: 28, - 3: 31, - 4: 30, - 5: 31, - 6: 30, - 7: 31, - 8: 31, - 9: 30, - 10: 31, - 11: 30, - 12: 31, - } - return daysinmonth_dict[month] - - def hypsometrystats(hyps_table, thickness_table): """Calculate the volume and mean associated with the hypsometry data. diff --git a/pygem/setup/config.py b/pygem/setup/config.py index c60f520c..b6bf47e4 100644 --- a/pygem/setup/config.py +++ b/pygem/setup/config.py @@ -131,7 +131,7 @@ def _validate_config(self, config): 'setup.rgi_region02': (str, type(None)), 'setup.glac_no_skip': (list, type(None)), 'setup.glac_no': (list, type(None)), - 'setup.min_glac_area_km2': int, + 'setup.min_glac_area_km2': (int, float), 'setup.include_landterm': bool, 'setup.include_laketerm': bool, 'setup.include_tidewater': bool, @@ -176,21 +176,21 @@ def _validate_config(self, config): 'calib.option_calibration': str, 'calib.priors_reg_fn': str, 'calib.HH2015_params': dict, - 'calib.HH2015_params.tbias_init': int, - 'calib.HH2015_params.tbias_step': int, - 'calib.HH2015_params.kp_init': float, - 'calib.HH2015_params.kp_bndlow': float, - 'calib.HH2015_params.kp_bndhigh': int, - 'calib.HH2015_params.ddfsnow_init': float, - 'calib.HH2015_params.ddfsnow_bndlow': float, - 'calib.HH2015_params.ddfsnow_bndhigh': float, + 'calib.HH2015_params.tbias_init': (int, float), + 'calib.HH2015_params.tbias_step': (int, float), + 'calib.HH2015_params.kp_init': (int, float), + 'calib.HH2015_params.kp_bndlow': (int, float), + 'calib.HH2015_params.kp_bndhigh': (int, float), + 'calib.HH2015_params.ddfsnow_init': (int, float), + 'calib.HH2015_params.ddfsnow_bndlow': (int, float), + 'calib.HH2015_params.ddfsnow_bndhigh': (int, float), 'calib.HH2015mod_params': dict, - 'calib.HH2015mod_params.tbias_init': int, - 'calib.HH2015mod_params.tbias_step': float, - 'calib.HH2015mod_params.kp_init': int, - 'calib.HH2015mod_params.kp_bndlow': float, - 'calib.HH2015mod_params.kp_bndhigh': int, - 'calib.HH2015mod_params.ddfsnow_init': float, + 'calib.HH2015mod_params.tbias_init': (int, float), + 'calib.HH2015mod_params.tbias_step': (int, float), + 'calib.HH2015mod_params.kp_init': (int, float), + 'calib.HH2015mod_params.kp_bndlow': (int, float), + 'calib.HH2015mod_params.kp_bndhigh': (int, float), + 'calib.HH2015mod_params.ddfsnow_init': (int, float), 'calib.HH2015mod_params.method_opt': str, 'calib.HH2015mod_params.params2opt': list, 'calib.HH2015mod_params.ftol_opt': float, @@ -199,22 +199,22 @@ def _validate_config(self, config): 'calib.emulator_params.emulator_sims': int, 'calib.emulator_params.overwrite_em_sims': bool, 'calib.emulator_params.opt_hh2015_mod': bool, - 'calib.emulator_params.tbias_step': float, - 'calib.emulator_params.tbias_init': int, - 'calib.emulator_params.kp_init': int, - 'calib.emulator_params.kp_bndlow': float, - 'calib.emulator_params.kp_bndhigh': int, - 'calib.emulator_params.ddfsnow_init': float, + 'calib.emulator_params.tbias_step': (int, float), + 'calib.emulator_params.tbias_init': (int, float), + 'calib.emulator_params.kp_init': (int, float), + 'calib.emulator_params.kp_bndlow': (int, float), + 'calib.emulator_params.kp_bndhigh': (int, float), + 'calib.emulator_params.ddfsnow_init': (int, float), 'calib.emulator_params.option_areaconstant': bool, 'calib.emulator_params.tbias_disttype': str, - 'calib.emulator_params.tbias_sigma': int, - 'calib.emulator_params.kp_gamma_alpha': int, - 'calib.emulator_params.kp_gamma_beta': int, + 'calib.emulator_params.tbias_sigma': (int, float), + 'calib.emulator_params.kp_gamma_alpha': (int, float), + 'calib.emulator_params.kp_gamma_beta': (int, float), 'calib.emulator_params.ddfsnow_disttype': str, - 'calib.emulator_params.ddfsnow_mu': float, - 'calib.emulator_params.ddfsnow_sigma': float, - 'calib.emulator_params.ddfsnow_bndlow': int, - 'calib.emulator_params.ddfsnow_bndhigh': float, + 'calib.emulator_params.ddfsnow_mu': (int, float), + 'calib.emulator_params.ddfsnow_sigma': (int, float), + 'calib.emulator_params.ddfsnow_bndlow': (int, float), + 'calib.emulator_params.ddfsnow_bndhigh': (int, float), 'calib.emulator_params.method_opt': str, 'calib.emulator_params.params2opt': list, 'calib.emulator_params.ftol_opt': float, @@ -222,33 +222,33 @@ def _validate_config(self, config): 'calib.MCMC_params': dict, 'calib.MCMC_params.option_use_emulator': bool, 'calib.MCMC_params.emulator_sims': int, - 'calib.MCMC_params.tbias_step': float, - 'calib.MCMC_params.tbias_stepsmall': float, + 'calib.MCMC_params.tbias_step': (int, float), + 'calib.MCMC_params.tbias_stepsmall': (int, float), 'calib.MCMC_params.option_areaconstant': bool, - 'calib.MCMC_params.mcmc_step': float, + 'calib.MCMC_params.mcmc_step': (int, float), 'calib.MCMC_params.n_chains': int, 'calib.MCMC_params.mcmc_sample_no': int, 'calib.MCMC_params.mcmc_burn_pct': int, 'calib.MCMC_params.thin_interval': int, 'calib.MCMC_params.ddfsnow_disttype': str, - 'calib.MCMC_params.ddfsnow_mu': float, - 'calib.MCMC_params.ddfsnow_sigma': float, - 'calib.MCMC_params.ddfsnow_bndlow': int, - 'calib.MCMC_params.ddfsnow_bndhigh': float, + 'calib.MCMC_params.ddfsnow_mu': (int, float), + 'calib.MCMC_params.ddfsnow_sigma': (int, float), + 'calib.MCMC_params.ddfsnow_bndlow': (int, float), + 'calib.MCMC_params.ddfsnow_bndhigh': (int, float), 'calib.MCMC_params.kp_disttype': str, 'calib.MCMC_params.tbias_disttype': str, - 'calib.MCMC_params.tbias_mu': int, - 'calib.MCMC_params.tbias_sigma': int, - 'calib.MCMC_params.tbias_bndlow': int, - 'calib.MCMC_params.tbias_bndhigh': int, - 'calib.MCMC_params.kp_gamma_alpha': int, - 'calib.MCMC_params.kp_gamma_beta': int, - 'calib.MCMC_params.kp_lognorm_mu': int, - 'calib.MCMC_params.kp_lognorm_tau': int, - 'calib.MCMC_params.kp_mu': int, - 'calib.MCMC_params.kp_sigma': float, - 'calib.MCMC_params.kp_bndlow': float, - 'calib.MCMC_params.kp_bndhigh': float, + 'calib.MCMC_params.tbias_mu': (int, float), + 'calib.MCMC_params.tbias_sigma': (int, float), + 'calib.MCMC_params.tbias_bndlow': (int, float), + 'calib.MCMC_params.tbias_bndhigh': (int, float), + 'calib.MCMC_params.kp_gamma_alpha': (int, float), + 'calib.MCMC_params.kp_gamma_beta': (int, float), + 'calib.MCMC_params.kp_lognorm_mu': (int, float), + 'calib.MCMC_params.kp_lognorm_tau': (int, float), + 'calib.MCMC_params.kp_mu': (int, float), + 'calib.MCMC_params.kp_sigma': (int, float), + 'calib.MCMC_params.kp_bndlow': (int, float), + 'calib.MCMC_params.kp_bndhigh': (int, float), 'calib.MCMC_params.option_calib_elev_change_1d': bool, 'calib.MCMC_params.rhoabl_disttype': str, 'calib.MCMC_params.rhoabl_mu': (int, float), @@ -285,26 +285,26 @@ def _validate_config(self, config): 'sim.out.export_extra_vars': bool, 'sim.out.export_binned_data': bool, 'sim.out.export_binned_components': bool, - 'sim.out.export_binned_area_threshold': int, + 'sim.out.export_binned_area_threshold': (int, float), 'sim.oggm_dynamics': dict, 'sim.oggm_dynamics.cfl_number': float, 'sim.oggm_dynamics.cfl_number_calving': float, 'sim.oggm_dynamics.glen_a_regional_relpath': str, 'sim.oggm_dynamics.use_regional_glen_a': bool, - 'sim.oggm_dynamics.fs': int, - 'sim.oggm_dynamics.glen_a_multiplier': int, - 'sim.icethickness_advancethreshold': int, - 'sim.terminus_percentage': int, + 'sim.oggm_dynamics.fs': (int, float), + 'sim.oggm_dynamics.glen_a_multiplier': (int, float), + 'sim.icethickness_advancethreshold': (int, float), + 'sim.terminus_percentage': (int, float), 'sim.params': dict, 'sim.params.use_constant_lapserate': bool, - 'sim.params.kp': int, - 'sim.params.tbias': int, - 'sim.params.ddfsnow': float, - 'sim.params.ddfsnow_iceratio': float, - 'sim.params.precgrad': float, - 'sim.params.lapserate': float, - 'sim.params.tsnow_threshold': int, - 'sim.params.calving_k': float, + 'sim.params.kp': (int, float), + 'sim.params.tbias': (int, float), + 'sim.params.ddfsnow': (int, float), + 'sim.params.ddfsnow_iceratio': (int, float), + 'sim.params.precgrad': (int, float), + 'sim.params.lapserate': (int, float), + 'sim.params.tsnow_threshold': (int, float), + 'sim.params.calving_k': (int, float), 'mb': dict, 'mb.option_surfacetype_initial': int, 'mb.include_firn': bool, @@ -321,12 +321,12 @@ def _validate_config(self, config): 'mb.Woodard_rf_opts.rf_month': int, 'mb.HH2015_rf_opts': dict, 'mb.HH2015_rf_opts.rf_layers': int, - 'mb.HH2015_rf_opts.rf_dz': int, - 'mb.HH2015_rf_opts.rf_dsc': int, - 'mb.HH2015_rf_opts.rf_meltcrit': float, - 'mb.HH2015_rf_opts.pp': float, - 'mb.HH2015_rf_opts.rf_dens_top': int, - 'mb.HH2015_rf_opts.rf_dens_bot': int, + 'mb.HH2015_rf_opts.rf_dz': (int, float), + 'mb.HH2015_rf_opts.rf_dsc': (int, float), + 'mb.HH2015_rf_opts.rf_meltcrit': (int, float), + 'mb.HH2015_rf_opts.pp': (int, float), + 'mb.HH2015_rf_opts.rf_dens_top': (int, float), + 'mb.HH2015_rf_opts.rf_dens_bot': (int, float), 'mb.HH2015_rf_opts.option_rf_limit_meltsnow': int, 'rgi': dict, 'rgi.rgi_relpath': str, @@ -346,13 +346,13 @@ def _validate_config(self, config): 'time.summer_month_start': int, 'time.timestep': str, 'constants': dict, - 'constants.density_ice': int, - 'constants.density_water': int, - 'constants.k_ice': float, - 'constants.k_air': float, - 'constants.ch_ice': int, - 'constants.ch_air': int, - 'constants.Lh_rf': int, + 'constants.density_ice': (int, float), + 'constants.density_water': (int, float), + 'constants.k_ice': (int, float), + 'constants.k_air': (int, float), + 'constants.ch_ice': (int, float), + 'constants.ch_air': (int, float), + 'constants.Lh_rf': (int, float), 'constants.tolerance': float, 'debug': dict, 'debug.refreeze': bool, diff --git a/pygem/setup/config.yaml b/pygem/setup/config.yaml index 83bf12b8..5697e790 100644 --- a/pygem/setup/config.yaml +++ b/pygem/setup/config.yaml @@ -212,13 +212,13 @@ calib: h_ref_relpath: /IceThickness_Farinotti/composite_thickness_RGI60-all_regions/ # 1d elevation change elev_change_1d: - elev_change_1d_relpath: /elev_change_1d/ # relative to main data path. per-glacier files within will be expected as _elev_change_1d_.json (e.g., 1.00570_elev_change_1d_.json) + elev_change_1d_relpath: /elev_change_1d/ # relative to main data path. per-glacier files within will be expected as _elev_change_1d_.json (e.g., 01.00570_elev_change_1d.json) # 1d melt extents meltextent_1d: - meltextent_1d_relpath: /SAR_data/merged/ # relative to main data path. per-glacier files within will be expected as _melt_extent_elev.csv (e.g., 01.00570_melt_extent_elev.csv) + meltextent_1d_relpath: /SAR_data/merged/ # relative to main data path. per-glacier files within will be expected as _melt_extent_elev.csv (e.g., 01.00570_melt_extent_elev.csv) # 1d snowlines snowline_1d: - snowline_1d_relpath: /SAR_data/merged/ # relative to main data path. per-glacier files within will be expected as _snowline_elev.csv (e.g., 01.00570_snowline_elev.csv) + snowline_1d_relpath: /SAR_data/merged/ # relative to main data path. per-glacier files within will be expected as _snowline_elev.csv (e.g., 01.00570_snowline_elev.csv) icethickness_cal_frac_byarea: 0.9 # Regional glacier area fraction that is used to calibrate the ice thickness # e.g., 0.9 means only the largest 90% of glaciers by area will be used to calibrate @@ -240,7 +240,7 @@ sim: - mad # export options (booleans) export_all_simiters: false # Exprort individual simulation results (false exports median and MAD from all sim_iters) - export_extra_vars: false # Option to export extra variables (temp, prec, melt, acc, etc.) + export_extra_vars: false # Option to export extra variables (temp, prec, melt, acc, etc.) export_binned_data: false # Export binned ice thickness export_binned_components: false # Export binned mass balance components (accumulation, melt, refreeze) export_binned_area_threshold: 0 # Area threshold for exporting binned ice thickness diff --git a/pygem/shop/mbdata.py b/pygem/shop/mbdata.py index 1ff56333..a507dfba 100755 --- a/pygem/shop/mbdata.py +++ b/pygem/shop/mbdata.py @@ -66,7 +66,6 @@ def mb_df_to_gdir( mbdata_fp = mbdata_fp + pygem_prms['calib']['data']['massbalance']['hugonnet2021_fn'] assert os.path.exists(mbdata_fp), 'Error, mass balance dataset does not exist: {mbdata_fp}' - assert 'hugonnet2021' in mbdata_fp.lower(), 'Error, mass balance dataset not yet supported: {mbdata_fp}' rgiid_cn = 'rgiid' mb_cn = 'mb_mwea' mberr_cn = 'mb_mwea_err' diff --git a/pygem/tests/test_04_auxiliary.py b/pygem/tests/test_04_auxiliary.py new file mode 100644 index 00000000..c9db58fc --- /dev/null +++ b/pygem/tests/test_04_auxiliary.py @@ -0,0 +1,87 @@ +import glob +import os +import subprocess + +import pytest + +from pygem.setup.config import ConfigManager + +""" +Test suite to any necessary aux""" + + +@pytest.fixture(scope='module') +def rootdir(): + config_manager = ConfigManager() + pygem_prms = config_manager.read_config() + return pygem_prms['root'] + + +def test_simulation_massredistribution_dynamics(rootdir): + """ + Test the run_simulation CLI script with the "MassRedistributionCurves" dynamical option. + """ + + # Run run_simulation CLI script + subprocess.run( + [ + 'run_simulation', + '-rgi_glac_number', + '1.03622', + '-option_calibration', + 'MCMC', + '-sim_climate_name', + 'ERA5', + '-sim_startyear', + '2000', + '-sim_endyear', + '2019', + '-nsims', + '1', + '-option_dynamics', + 'MassRedistributionCurves', + '-outputfn_sfix', + 'mrcdynamics_', + ], + check=True, + ) + + # Check if output files were created + outdir = os.path.join(rootdir, 'Output', 'simulations', '01', 'ERA5') + output_files = glob.glob(os.path.join(outdir, '**', '*_mrcdynamics_all.nc'), recursive=True) + assert output_files, f'Simulation output file not found in {outdir}' + + +def test_simulation_no_dynamics(rootdir): + """ + Test the run_simulation CLI script with no dynamics option. + """ + + # Run run_simulation CLI script + subprocess.run( + [ + 'run_simulation', + '-rgi_glac_number', + '1.03622', + '-option_calibration', + 'MCMC', + '-sim_climate_name', + 'ERA5', + '-sim_startyear', + '2000', + '-sim_endyear', + '2019', + '-nsims', + '1', + '-option_dynamics', + 'None', + '-outputfn_sfix', + 'nodynamics_', + ], + check=True, + ) + + # Check if output files were created + outdir = os.path.join(rootdir, 'Output', 'simulations', '01', 'ERA5') + output_files = glob.glob(os.path.join(outdir, '**', '*_nodynamics_all.nc'), recursive=True) + assert output_files, f'Simulation output file not found in {outdir}' diff --git a/pygem/tests/test_04_postproc.py b/pygem/tests/test_05_postproc.py similarity index 78% rename from pygem/tests/test_04_postproc.py rename to pygem/tests/test_05_postproc.py index 3d50aa02..35a60eb8 100644 --- a/pygem/tests/test_04_postproc.py +++ b/pygem/tests/test_05_postproc.py @@ -16,14 +16,24 @@ def rootdir(): return pygem_prms['root'] -def test_postproc_monthly_mass(rootdir): +def test_postproc_subannual_mass(rootdir): """ - Test the postproc_monthly_mass CLI script. + Test the postproc_subannual_mass CLI script. """ simdir = os.path.join(rootdir, 'Output', 'simulations', '01', 'CESM2', 'ssp245', 'stats') # Run postproc_monthyl_mass CLI script - subprocess.run(['postproc_monthly_mass', '-simdir', simdir], check=True) + subprocess.run(['postproc_subannual_mass', '-simdir', simdir], check=True) + + +def test_postproc_binned_subannual_thick(rootdir): + """ + Test the postproc_binned_subannual_thick CLI script. + """ + simdir = os.path.join(rootdir, 'Output', 'simulations', '01', 'CESM2', 'ssp245', 'binned') + + # Run postproc_monthyl_mass CLI script + subprocess.run(['postproc_binned_subannual_thick', '-simdir', simdir], check=True) def test_postproc_compile_simulations(rootdir): @@ -63,12 +73,12 @@ def test_check_compiled_product(rootdir): """ # skip variables that are not in the compiled products vars_to_skip = [ - 'glac_temp_monthly', + 'glac_temp', 'glac_mass_change_ignored_annual', - 'offglac_prec_monthly', - 'offglac_refreeze_monthly', - 'offglac_melt_monthly', - 'offglac_snowpack_monthly', + 'offglac_prec', + 'offglac_refreeze', + 'offglac_melt', + 'offglac_snowpack', ] simpath = os.path.join( @@ -79,7 +89,7 @@ def test_check_compiled_product(rootdir): 'CESM2', 'ssp245', 'stats', - '1.03622_CESM2_ssp245_MCMC_ba1_50sets_2000_2100_all.nc', + '1.03622_CESM2_ssp245_MCMC_ba1_10sets_2000_2100_all.nc', ) compdir = os.path.join(rootdir, 'Output', 'simulations', 'compile', 'glacier_stats') @@ -104,13 +114,13 @@ def test_check_compiled_product(rootdir): # pull data values simvals = simvar.values - compvals = compvar.values[0, :, :] # first index is the glacier index + compvals = compvar.values[0, :, :] # 0th index is the glacier index # check that compiled product has same shape as original data assert simvals.shape == compvals.shape, ( f'Compiled product shape {compvals.shape} does not match original data shape {simvals.shape}' ) # check that compiled product matches original data - assert np.all(np.array_equal(simvals, compvals)), ( + assert np.allclose(simvals, compvals, rtol=1e-8, atol=1e-12, equal_nan=True), ( f'Compiled product for {var} does not match original data' ) diff --git a/pygem/utils/_funcs.py b/pygem/utils/_funcs.py index 84ff15eb..ae57bf8a 100755 --- a/pygem/utils/_funcs.py +++ b/pygem/utils/_funcs.py @@ -52,14 +52,14 @@ def annualweightedmean_array(var, dates_table): var : np.ndarray Variable with monthly or daily timestep dates_table : pd.DataFrame - Table of dates, year, month, daysinmonth, wateryear, and season for each timestep + Table of dates, year, month, days_in_step, wateryear, and season for each timestep Returns ------- var_annual : np.ndarray Annual weighted mean of variable """ if pygem_prms['time']['timestep'] == 'monthly': - dayspermonth = dates_table['daysinmonth'].values.reshape(-1, 12) + dayspermonth = dates_table['days_in_step'].values.reshape(-1, 12) # creates matrix (rows-years, columns-months) of the number of days per month daysperyear = dayspermonth.sum(axis=1) # creates an array of the days per year (includes leap years) @@ -76,10 +76,10 @@ def annualweightedmean_array(var, dates_table): if var_annual.shape[1] == 1: var_annual = var_annual.reshape(var_annual.shape[0]) elif pygem_prms['time']['timestep'] == 'daily': - print( - '\nError: need to code the groupbyyearsum and groupbyyearmean for daily timestep.Exiting the model run.\n' - ) - exit() + var_annual = var.mean(1) + else: + # var_annual = var.mean(1) + assert 1 == 0, 'add this functionality for weighting that is not monthly or daily' return var_annual diff --git a/pyproject.toml b/pyproject.toml index 7ce3a44e..4d6b03ba 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -50,8 +50,8 @@ run_calibration_frontalablation = "pygem.bin.run.run_calibration_frontalablation run_calibration = "pygem.bin.run.run_calibration:main" run_mcmc_priors = "pygem.bin.run.run_mcmc_priors:main" run_simulation = "pygem.bin.run.run_simulation:main" -postproc_monthly_mass = "pygem.bin.postproc.postproc_monthly_mass:main" -postproc_binned_monthly_mass = "pygem.bin.postproc.postproc_binned_monthly_mass:main" +postproc_subannual_mass = "pygem.bin.postproc.postproc_subannual_mass:main" +postproc_binned_subannual_thick = "pygem.bin.postproc.postproc_binned_subannual_thick:main" postproc_distribute_ice = "pygem.bin.postproc.postproc_distribute_ice:main" postproc_compile_simulations = "pygem.bin.postproc.postproc_compile_simulations:main" list_failed_simulations = "pygem.bin.op.list_failed_simulations:main"