diff --git a/docs/install_pygem.md b/docs/install_pygem.md index f645c9d8..87c92f33 100644 --- a/docs/install_pygem.md +++ b/docs/install_pygem.md @@ -6,7 +6,6 @@ The model is stored in two repositories [(see model structure)](model_structure_ A conda environment is a directory that contains a specific collection of installed packages. The use of environments reduces issues caused by package dependencies. The model is designed to be compatible with OGGM. We therefore get started by following the [installation instructions from OGGM](https://docs.oggm.org/en/stable/installing-oggm.html). Once your conda environment is setup for OGGM, add the core of PyGEM using pip. - ``` pip install pygem ``` @@ -14,14 +13,16 @@ pip install pygem This will provide you with a conda environment that has the basic functionality to run simple calibration options (e.g., 'HH2015', 'HH2015mod') and simulations. If you want to use the emulators and Bayesian inference, the advanced environment is required. ### Developing PyGEM -Are you interested in developing PyGEM? If so, we recommend uninstalling pygem, forking the [PyGEM's github repository](https://github.com/drounce/PyGEM) and then [cloning the github repository](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) onto your local machine. +Are you interested in developing PyGEM? If so, we recommend forking the [PyGEM's github repository](https://github.com/PyGEM-Community/PyGEM) and then [cloning the github repository](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) onto your local machine. +Note, if PyGEM was already installed via PyPI, first uninstall: ``` pip uninstall pygem -git clone https://github.com/drounce/PyGEM.git +```` + +You can then use pip to install your locally cloned fork of PyGEM in 'editable' mode like so: ``` -```{warning} -The directory structure here is important. The cloned PyGEM directory should be on the same level as the PyGEM-scripts. +pip install -e /path/to/your/PyGEM/clone ``` ### Advanced environment: GPyTorch (emulator only) @@ -64,7 +65,7 @@ The only way to find out if your package dependencies work is to test it by runn ## Install PyGEM-Scripts -The scripts that are used to run PyGEM are located in the [PyGEM-Scripts repository](https://github.com/drounce/PyGEM-scripts) on github. To run the model, you can either (i) clone the repository or (ii) fork the repository to develop/add your own scripts. For instructions, follow github’s instructions on [cloning](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) or [forking a repository](https://docs.github.com/en/get-started/quickstart/fork-a-repo). Once the repository is installed on your local machine, you can run the model from this directory. +The scripts that are used to run PyGEM are located in the [PyGEM-Scripts repository](https://github.com/PyGEM-Community/PyGEM-scripts) on github. To run the model, you can either (i) clone the repository or (ii) fork the repository to develop/add your own scripts. For instructions, follow github’s instructions on [cloning](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) or [forking a repository](https://docs.github.com/en/get-started/quickstart/fork-a-repo). Once the repository is installed on your local machine, you can run the model from this directory. ```{note} Be sure that your [directory structure](directory_structure_target) is setup properly before you try running the model! diff --git a/pygem/gcmbiasadj.py b/pygem/gcmbiasadj.py index 6419e54e..2a2d6892 100755 --- a/pygem/gcmbiasadj.py +++ b/pygem/gcmbiasadj.py @@ -321,8 +321,9 @@ def prec_biasadj_opt1(ref_prec, ref_elev, gcm_prec, dates_table_ref, dates_table gcm_prec_biasadj_subset = ( gcm_prec_biasadj[:,gcm_subset_idx_start:gcm_subset_idx_end+1][:,gcm_spinupyears*12:]) gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / ref_prec_nospinup.sum(axis=1) - assert np.min(gcm_prec_biasadj_frac) > 0.5 and np.max(gcm_prec_biasadj_frac) < 2, ( - 'Error with gcm precipitation bias adjustment: total ref and gcm prec differ by more than factor of 2') + assert np.min(gcm_prec_biasadj_frac) > 0.45 and np.max(gcm_prec_biasadj_frac) < 2, ( + f'Error with gcm precipitation bias adjustment: total ref and gcm prec differ by more than factor of 2.' + f' min: {np.min(gcm_prec_biasadj_frac)}, max: {np.max(gcm_prec_biasadj_frac)}') assert gcm_prec_biasadj.max() <= 10, 'gcm_prec_adj (precipitation bias adjustment) too high, needs to be modified' assert gcm_prec_biasadj.min() >= 0, 'gcm_prec_adj is producing a negative precipitation value' diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index fc6d7324..c79132e1 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -21,6 +21,303 @@ cfg.initialize() +# calving law + #%% ----- Calving Law ---------- + #%% ----- Calving Law ---------- +def fa_sermeq_speed_law(model,last_above_wl, v_scaling=1, verbose=False, + tau0=1.5, yield_type='constant', mu=0.01, + trim_profile=0): + """ + This function is used to calculate frontal ablation given ice speed forcing, + for lake-terminating and tidewater glaciers + + @author: Ruitang Yang & Lizz Ultee + + Authors: Ruitang Yang & Lizz Ultee + Parameters + ---------- + + model : oggm.core.flowline.FlowlineModel + the model instance calling the function + flowline : oggm.core.flowline.Flowline + the instance of the flowline object on which the calving law is called + fl_id : float, optional + the index of the flowline in the fls array (might be ignored by some MB models) + last_above_wl : int + the index of the last pixel above water (in case you need to know + where it is). + v_scaling: float + velocity scaling factor, >0, default is 1 + Terminus_mb : array + Mass balance along the flowline or nearest the terminus [m/a]. Default None... + TODO: set default behavior, check the unit meter of ice per year or m w.e. per year? + verbose: Boolean, optional + Whether to print component parts for inspection. Default False. + + tau0: float, optional + This glacier's yield strength [Pa]. Default is 150 kPa. + yield_type: str, optional + 'constant' or 'variable' (Mohr-Coulomb) yielding. Default is constant. + mu: float, optional + Mohr-Coulomb cohesion, a coefficient between 0 and 1. Default is 0.01. + Only used if we have variable yield + + trim_profile: int, optional + How many grid cells at the end of the profile to ignore. Default is 1. + If the initial profile is set by k-calving (as in testing) there can be a + weird cliff shape with very thin final grid point and large velocity gradient + + Returns + ------- + fa_viscoplastic: float + Frontal ablation rate [m/a] based on viscoplastic assumptions + SQFA: dict + Frontal ablation rate [m/a] based on viscoplastic assumptions + serface elevation at the terminus [m a.s.l.] based on OGGM ## TODO CHANGE IT BY PyGEM Surface mass balance result + bed elevation at the terminus [m a.s.l.] + Terminus Thickness [m] + Yield terminus thickness [m] + Velocity at the terminus [m/a] + Surface mass balance at the terminus [m/a] + Length change at the terminus [m/a] based on viscoplastic assumptions + TODO: output the length change in case we + have the length change results from observations in-situ or remote sensing (Thomas Schellenberger has the machine learning products) + Frontal ablation rate [m/a] based on viscoplastic assumptions + + Explanation of sign + ------- + fa_viscoplastic: negative, mass loss + dLdt: length change rate, positive if advance; negative if retreat + terminus mass balance: negative if mass loss; positive if mass gain + """ + # --------------------------------------------------------------------------- + # class NegativeValueError(Exception): + # pass + # --------------------------------------------------------------------------- + ## Global constants + G = 9.8 # acceleration due to gravity in m/s^2 + RHO_ICE = 920.0 # ice density kg/m^3 + RHO_SEA = 1020.0 # seawater density kg/m^3 + + # --------------------------------------------------------------------------- + # the yield strength + def tau_y(tau0=1.5, yield_type='constant', bed_elev=None, thick=None, mu=0.01): + """ + Functional form of yield strength. + Can do constant or Mohr-Coulomb yield strength. Ideally, the glacier's yield type + ('constant' or 'variable') would be saved in a model instance. + + Parameters + ---------- + tau0: float, optional + Initial guess for yield strength [Pa]. Default is 150 kPa. + yield_type: str, optional + 'constant' or 'variable' (Mohr-Coulomb) yielding. Default is constant. + bed_elev: float, optional + Bed elevation, dimensional [m]. The default is None. + thick: float, optional + Ice thickness, dimensional [m]. The default is None. + mu: float, optional + Mohr-Coulomb cohesion, a coefficient between 0 and 1. Default is 0.01. + + Returns + ------- + tau_y: float + The yield strength for these conditions. + """ + tau1=tau0*1e5 + if yield_type == 'variable': + try: + if bed_elev < 0: + D = -1 * bed_elev # Water depth D the nondim bed topography value when Z<0 + else: + D = 0 + except: + print('You must set a bed elevation and ice thickness to use variable yield strength.') + N = RHO_ICE * G * thick - RHO_SEA * G * D # Normal stress at bed + ty = tau1 + mu * N + else: # assume constant if not set + ty = tau1 + return ty + + + # --------------------------------------------------------------------------- + # calculate the yield ice thickness + + def balance_thickness(yield_strength, bed_elev): + """ + Ice thickness such that the stress matches the yield strength. + + Parameters + ---------- + yield_strength: float + The yield strength near the terminus. + If yield type is constant, this will of course be the same everywhere. If yield type is + variable (Mohr-Coulomb), the yield strength at the terminus could differ from elsewhere. + bed_elev: float + Elevation of glacier bed at the terminus + + Returns + ------- + Hy: float + The ice thickness for stress balance at the terminus. + """ + if bed_elev < 0: + D = -1 * bed_elev + else: + D = 0 + return (2 * yield_strength / (RHO_ICE * G)) + np.sqrt( + (RHO_SEA * (D ** 2) / RHO_ICE) + ((2 * yield_strength / (RHO_ICE * G)) ** 2)) + # TODO: Check on exponent on last term. In Ultee & Bassis 2016, this is squared, but in Ultee & Bassis 2020 supplement, it isn't. + + # --------------------------------------------------------------------------- + # calculate frontal ablation based on the ice thickness, speed at the terminus + fls=model.fls + flowline=fls[-1] + surface_m = flowline.surface_h + bed_m = flowline.bed_h + width_m = flowline.widths_m + velocity_m = model.u_stag[-1]*cfg.SEC_IN_YEAR + x_m = flowline.dis_on_line*flowline.map_dx/1000 + + # gdir : py:class:`oggm.GlacierDirectory` + # the glacier directory to process + # fls = model.gdir.read_pickle('model_flowlines') + # mbmod_fl = massbalance.MultipleFlowlineMassBalance(model.gdir, fls=fls, use_inversion_flowlines=True, + # mb_model_class=MonthlyTIModel) + mb_annual=model.mb_model.get_annual_mb(heights=surface_m, fl_id=-1, year=model.yr, fls=model.fls) + + Terminus_mb = mb_annual*cfg.SEC_IN_YEAR + # slice up to index+1 to include the last nonzero value + # profile: NDarray + # The current profile (x, surface, bed,width) as calculated by the base model + # Unlike core SERMeQ, these should be DIMENSIONAL [m]. + profile=(x_m[:last_above_wl+1], + surface_m[:last_above_wl+1], + bed_m[:last_above_wl+1],width_m[:last_above_wl+1]) + # model_velocity: array + # Velocity along the flowline [m/a] as calculated by the base model + # Should have values for the points nearest the terminus...otherwise + # doesn't matter if this is the same shape as the profile array. + # TODO: Check with the remote sensing products, or at least to validate the model products + model_velocity=velocity_m[:last_above_wl+1] + # remove lowest cells if needed + last_index = -1 * (trim_profile + 1) + ## TODO: Check the flowline model, the decrease the distance between two adjacent points along the flowline, and then calculate the averaged gradient for dhdx,dhydx,dudx + ## + if isinstance(Terminus_mb, (int, float)): + terminus_mb = Terminus_mb + elif isinstance(Terminus_mb, (list, np.ndarray)): + terminus_mb = Terminus_mb[last_index] + else: + print("please input the correct mass balance datatype") + # + if isinstance(model_velocity, (int, float)): + model_velocity = v_scaling * model_velocity + elif isinstance(model_velocity, list): + model_velocity = v_scaling * np.array(model_velocity) + elif isinstance(model_velocity, np.ndarray): + model_velocity = v_scaling * model_velocity + else: + print("please input the correct velocity datatype") + ## Ice thickness and yield thickness nearest the terminus + se_terminus = profile[1][last_index] + bed_terminus = profile[2][last_index] + h_terminus = se_terminus - bed_terminus + width_terminus = profile[3][last_index] + tau_y_terminus = tau_y(tau0=tau0, bed_elev=bed_terminus, thick=h_terminus, yield_type=yield_type) + Hy_terminus = balance_thickness(yield_strength=tau_y_terminus, bed_elev=bed_terminus) + if isinstance(model_velocity, (int, float)): + U_terminus = model_velocity + U_adj = model_velocity + else: + U_terminus = model_velocity[last_index] ## velocity, assuming last point is terminus + U_adj = model_velocity[last_index - 1] + ## Ice thickness and yield thickness at adjacent point + se_adj = profile[1][last_index - 1] + bed_adj = profile[2][last_index - 1] + H_adj = se_adj - bed_adj + tau_y_adj = tau_y(tau0=tau0, bed_elev=bed_adj, thick=H_adj, yield_type=yield_type) + Hy_adj = balance_thickness(yield_strength=tau_y_adj, bed_elev=bed_adj) + # Gradients + dx_term = profile[0][last_index] - profile[0][last_index - 1] ## check grid spacing close to terminus + dHdx = (h_terminus - H_adj) / dx_term + dHydx = (Hy_terminus - Hy_adj) / dx_term + if np.isnan(U_terminus) or np.isnan(U_adj): + dUdx = np.nan ## velocity gradient + ## Group the terms + dLdt_numerator = np.nan + dLdt_denominator = np.nan ## TODO: compute dHydx + dLdt_viscoplastic = np.nan + # fa_viscoplastic = dLdt_viscoplastic -U_terminus ## frontal ablation rate + fa_viscoplastic = np.nan ## frontal ablation rate + else: + # Gradients + # dx_term = profile[0][last_index] - profile[0][last_index - 1] ## check grid spacing close to terminus + # dHdx = (h_terminus - H_adj) / dx_term + # dHydx = (Hy_terminus - Hy_adj) / dx_term + dUdx = (U_terminus - U_adj) / dx_term ## velocity gradient + ## Group the terms + dLdt_numerator = terminus_mb - (h_terminus * dUdx) - (U_terminus * dHdx) + dLdt_denominator = dHydx - dHdx ## TODO: compute dHydx + dLdt_viscoplastic = dLdt_numerator / dLdt_denominator + # fa_viscoplastic = dLdt_viscoplastic -U_terminus ## frontal ablation rate + + # try: + U_calving = U_terminus - dLdt_viscoplastic ## frontal ablation rate + fa_viscoplastic=U_calving + # if U_calving<0: + # print("The glacier is advancing, and the advancing rate is larger than ice flow speed at the terminus, please check ") + # if U_calving>0 or U_calving==0: + # fa_viscoplastic=U_calving + # else: + # fa_viscoplastic=U_calving + # # fa_viscoplastic=np.nan + # raise NegativeValueError("Something is wrong, right now the calving in negative, which should be positive or zero") + # except NegativeValueError as e: + # print ("The glacier is advancing, and the advancing rate is larger than ice flow speed at the terminus, please check ") + + + SQFA = {'se_terminus': se_terminus, + 'bed_terminus': bed_terminus, + 'Thickness_termi': h_terminus, + 'Width_termi': width_terminus, + 'Hy_thickness': Hy_terminus, + 'Velocity_termi': U_terminus, + 'Terminus_mb': terminus_mb, + 'dLdt': dLdt_viscoplastic, + 'Sermeq_fa': fa_viscoplastic} + if verbose: + print('For inspection on debugging - all should be DIMENSIONAL (m/a):') + # print('profile_length={}'.format(profile_length)) + print('last_index={}'.format(last_index)) + print('se_terminus={}'.format(se_terminus)) + print('bed_terminus={}'.format(bed_terminus)) + print('se_adj={}'.format(se_adj)) + print('bed_adj={}'.format(bed_adj)) + print('Thicknesses: Hterm {}, Hadj {}'.format(h_terminus, H_adj)) + print('Hy_terminus={}'.format(Hy_terminus)) + print('Hy_adj={}'.format(Hy_adj)) + print('U_terminus={}'.format(U_terminus)) + print('U_adj={}'.format(U_adj)) + print('dUdx={}'.format(dUdx)) + print('dx_term={}'.format(dx_term)) + print('Checking dLdt: terminus_mb = {}. \n H dUdx = {}. \n U dHdx = {}.'.format(terminus_mb, dUdx * h_terminus, + U_terminus * dHdx)) + print('Denom: dHydx = {} \n dHdx = {}'.format(dHydx, dHdx)) + print('Viscoplastic dLdt={}'.format(dLdt_viscoplastic)) + print('Terminus surface mass balance ma= {}'.format(terminus_mb)) + print('Sermeq frontal ablation ma={}'.format(fa_viscoplastic)) + else: + pass + return SQFA + + + + + + + #%% class MassRedistributionCurveModel(FlowlineModel): @@ -455,8 +752,11 @@ def updategeometry(self, year, debug=False): 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() - - + + + + + #%% ----- FRONTAL ABLATION ----- def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, calving_k=None, debug=False ): @@ -520,7 +820,13 @@ def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, h = fl.thick[last_above_wl] d = h - (fl.surface_h[last_above_wl] - self.water_level) k = self.calving_k - q_calving = k * d * h * fl.widths_m[last_above_wl] + #q_calving = k * d * h * fl.widths_m[last_above_wl] + # Sermeq calving law + print('****************Sermeq claving law start********************') + s_fa = fa_sermeq_speed_law(self,last_above_wl=last_above_wl,v_scaling = 1, verbose = False,tau0 = self.calving_k, + yield_type = 'constant', mu = 0.01,trim_profile = 0) + q_calving = s_fa ['Sermeq_fa']*s_fa['Thickness_termi']*s_fa['Width_termi']/cfg.SEC_IN_YEAR + print('****************Sermeq claving law end successfully********************') # Max frontal ablation is removing all bins with bed below water level glac_idx_bsl = np.where((fl.thick > 0) & (fl.bed_h < self.water_level))[0] diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 3cb9347a..4badbbb5 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -7,6 +7,7 @@ """ # External libraries import numpy as np +import traceback # Local libraries from oggm.core.massbalance import MassBalanceModel import pygem_input as pygem_prms @@ -48,6 +49,10 @@ def __init__(self, gdir, modelprms, glacier_rgi_table, hindcast : Boolean switch to run the model in reverse or not (may be irrelevant after converting to OGGM's setup) """ + + #print("the input flowlines in PyGEMMassBalance Model is: ") + #print(type(fls).__name__) + if debug: print('\n\nDEBUGGING MASS BALANCE FUNCTION\n\n') self.debug_refreeze = debug_refreeze @@ -120,6 +125,7 @@ def __init__(self, gdir, modelprms, glacier_rgi_table, 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_lengthchange = np.zeros((nbins,self.nmonths)) self.glac_bin_snowpack = np.zeros((nbins,self.nmonths)) self.glac_bin_massbalclim = np.zeros((nbins,self.nmonths)) self.glac_bin_massbalclim_annual = np.zeros((nbins,self.nyears)) @@ -138,7 +144,9 @@ def __init__(self, gdir, modelprms, glacier_rgi_table, 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_length_change = np.zeros(self.nmonths) self.glac_wide_massbaltotal = np.zeros(self.nmonths) + self.glac_wide_massbalclim = np.zeros(self.nmonths) self.glac_wide_runoff = np.zeros(self.nmonths) self.glac_wide_snowline = np.zeros(self.nmonths) self.glac_wide_area_annual = np.zeros(self.nyears+1) @@ -180,9 +188,20 @@ def __init__(self, gdir, modelprms, glacier_rgi_table, self.sea_level = 0 rgi_region = int(glacier_rgi_table.RGIId.split('-')[1].split('.')[0]) - - def get_annual_mb(self, heights, year=None, fls=None, fl_id=None, + def get_monthly_mb(self, heights, year=None, fls=None, fl_id=None, debug=False, option_areaconstant=False): + print("year in get_monthly_mb is :",year) + year_floor=np.floor(year) + month=year-int(year_floor) + print("year_floor to get_annual_mb is :",year_floor) + print("month to get_annual_mb is :",month) + mb=self.get_annual_mb(heights=heights, year=year_floor, fls=fls, fl_id=fl_id, + debug=debug, option_areaconstant=False, year_month=month) + + return mb + + def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, + debug=False, option_areaconstant=False, year_month=None): """FIXED FORMAT FOR THE FLOWLINE MODEL Returns annual climatic mass balance [m ice per second] @@ -198,469 +217,522 @@ def get_annual_mb(self, heights, year=None, fls=None, fl_id=None, ------- mb : np.array mass balance for each bin [m ice per second] - """ - year = int(year) - if self.repeat_period: - year = year % (pygem_prms.gcm_endyear - pygem_prms.gcm_startyear) - - fl = fls[fl_id] - np.testing.assert_allclose(heights, fl.surface_h) - glacier_area_t0 = fl.widths_m * fl.dx_meter - glacier_area_initial = self.glacier_area_initial - fl_widths_m = getattr(fl, 'widths_m', None) - fl_section = getattr(fl,'section',None) - # Ice thickness (average) - if fl_section is not None and fl_widths_m is not None: - icethickness_t0 = np.zeros(fl_section.shape) - icethickness_t0[fl_widths_m > 0] = fl_section[fl_widths_m > 0] / fl_widths_m[fl_widths_m > 0] - else: - icethickness_t0 = None - - # Quality control: ensure you only have glacier area where there is ice - if icethickness_t0 is not None: - glacier_area_t0[icethickness_t0 == 0] = 0 + """ + #year = int(year) + year=int(year) + #print("Year is: ", year, "Month is:", year_month) + try: + if self.repeat_period: + year = year % (pygem_prms.gcm_endyear - pygem_prms.gcm_startyear) + # print("######################### In PyGEMMassBalance get annual mb") + # print("fl_id is ",fl_id) + # print("fls is ",fls) + # print("#########################") + # try: + # fls = self.gdir.read_pickle('inversion_flowlines') + # print("######################### 2") + # print("fl_id is ",fl_id) + # print("fls is ",fls) + # print("######################### 2") + # except: + # print("**************Something is wrong with the fls read**************") + # print(traceback.format_exc()) + + if (fls is not None) and (fl_id is not None): + #print("fls and fl_id are not None") + try: + fl = fls[fl_id] + except: + print(traceback.format_exc()) + else: + print("Please check the fls and fl_id, they should not be None") + print(traceback.format_exc()) + err_heights = heights - fl.surface_h + if max(abs(err_heights)) > 0: + print('WARNING: Heights do not match flowline surface heights') + print('Heights:', heights) + print('Flowline surface heights:', fl.surface_h) + print('Max absolute difference:', max(abs(err_heights))) + #raise ValueError('Heights do not match flowline surface heights') + np.testing.assert_allclose(heights, fl.surface_h) + glacier_area_t0 = fl.widths_m * fl.dx_meter + glacier_area_initial = self.glacier_area_initial + fl_widths_m = getattr(fl, 'widths_m', None) + fl_section = getattr(fl,'section',None) + # Ice thickness (average) + if fl_section is not None and fl_widths_m is not None: + icethickness_t0 = np.zeros(fl_section.shape) + icethickness_t0[fl_widths_m > 0] = fl_section[fl_widths_m > 0] / fl_widths_m[fl_widths_m > 0] + else: + icethickness_t0 = None + #print('******************** attributes of the flowline at the time t0 ******************** ') + #print('glacier_area_initial is:',glacier_area_initial) + #print('glacier_area_t0 is:',glacier_area_t0) + #print('the fl_widths_m at t0 is:',fl_widths_m) + #print('the fl_section area at t0 is:',fl_section) + #print('the ice thickness at t0 is:',icethickness_t0) + # Quality control: ensure you only have glacier area where there is ice + if icethickness_t0 is not None: + glacier_area_t0[icethickness_t0 == 0] = 0 + + # Record ice thickness + self.glac_bin_icethickness_annual[:,year] = icethickness_t0 - # Record ice thickness - self.glac_bin_icethickness_annual[:,year] = icethickness_t0 - - # Glacier indices - glac_idx_t0 = glacier_area_t0.nonzero()[0] - - nbins = heights.shape[0] - nmonths = self.glacier_gcm_temp.shape[0] - - # Local variables - bin_precsnow = np.zeros((nbins,nmonths)) - - # Refreezing specific layers - if pygem_prms.option_refreezing == 'HH2015' and year == 0: - self.te_rf[:,:,0] = 0 # layer temp of each elev bin for present time step - self.tl_rf[:,:,0] = 0 # layer temp of each elev bin for previous time step - elif pygem_prms.option_refreezing == 'Woodward': - 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 == 0: - self.surfacetype, self.firnline_idx = self._surfacetypebinsinitial(self.heights) - self.glac_bin_surfacetype_annual[:,year] = self.surfacetype - - # Off-glacier area and indices - if option_areaconstant == False: - self.offglac_bin_area_annual[:,year] = glacier_area_initial - glacier_area_t0 - offglac_idx = np.where(self.offglac_bin_area_annual[:,year] > 0)[0] - - # Functions currently set up for monthly timestep - # only compute mass balance while glacier exists - if (pygem_prms.timestep == 'monthly'): -# if (pygem_prms.timestep == 'monthly') and (glac_idx_t0.shape[0] != 0): - - # AIR TEMPERATURE: Downscale the gcm temperature [deg C] to each bin - if pygem_prms.option_temp2bins == 1: - # 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[:,12*year:12*(year+1)] = (self.glacier_gcm_temp[12*year:12*(year+1)] + - self.glacier_gcm_lrgcm[12*year:12*(year+1)] * - (self.glacier_rgi_table.loc[pygem_prms.option_elev_ref_downscale] - self.glacier_gcm_elev) + - self.glacier_gcm_lrglac[12*year:12*(year+1)] * (heights - - self.glacier_rgi_table.loc[pygem_prms.option_elev_ref_downscale])[:, np.newaxis] + - self.modelprms['tbias']) - - # PRECIPITATION/ACCUMULATION: Downscale the precipitation (liquid and solid) to each bin - if pygem_prms.option_prec2bins == 1: - # Precipitation using precipitation factor and precipitation gradient - # P_bin = P_gcm * prec_factor * (1 + prec_grad * (z_bin - z_ref)) - bin_precsnow[:,12*year:12*(year+1)] = (self.glacier_gcm_prec[12*year:12*(year+1)] * - self.modelprms['kp'] * (1 + self.modelprms['precgrad'] * (heights - - self.glacier_rgi_table.loc[pygem_prms.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.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,12*year:12*(year+1)] = ( - bin_precsnow[glac_idx_75,12*year:12*(year+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 - 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.option_accumulation == 1: - # if temperature above threshold, then rain - (self.bin_prec[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold']]) = ( - bin_precsnow[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold']]) - # if temperature below threshold, then snow - (self.bin_acc[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] <= self.modelprms['tsnow_threshold']]) = ( - bin_precsnow[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] <= self.modelprms['tsnow_threshold']]) - elif pygem_prms.option_accumulation == 2: - # if temperature between min/max, then mix of snow/rain using linear relationship between min/max - self.bin_prec[:,12*year:12*(year+1)] = ( - (0.5 + (self.bin_temp[:,12*year:12*(year+1)] - - self.modelprms['tsnow_threshold']) / 2) * bin_precsnow[:,12*year:12*(year+1)]) - self.bin_acc[:,12*year:12*(year+1)] = ( - bin_precsnow[:,12*year:12*(year+1)] - self.bin_prec[:,12*year:12*(year+1)]) - # if temperature above maximum threshold, then all rain - (self.bin_prec[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold'] + 1]) = ( - bin_precsnow[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold'] + 1]) - (self.bin_acc[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold'] + 1]) = 0 - # if temperature below minimum threshold, then all snow - (self.bin_acc[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] <= self.modelprms['tsnow_threshold'] - 1]) = ( - bin_precsnow[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] <= self.modelprms['tsnow_threshold'] - 1]) - (self.bin_prec[:,12*year:12*(year+1)] - [self.bin_temp[:,12*year:12*(year+1)] <= 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 = 12*year + 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.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.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]) - # 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.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.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.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.rf_dsc - - if pygem_prms.option_rf_limit_meltsnow == 1: - bin_meltlimit = self.bin_meltsnow.copy() + # Glacier indices + glac_idx_t0 = glacier_area_t0.nonzero()[0] + #print('the indices of glacier nonzero at t0 is:',glac_idx_t0) + + nbins = heights.shape[0] + nmonths = self.glacier_gcm_temp.shape[0] + #print('nbins is',nbins,'nmonths is :',nmonths) + # Local variables + bin_precsnow = np.zeros((nbins,nmonths)) + + # Refreezing specific layers + if pygem_prms.option_refreezing == 'HH2015' and year == 0: + self.te_rf[:,:,0] = 0 # layer temp of each elev bin for present time step + self.tl_rf[:,:,0] = 0 # layer temp of each elev bin for previous time step + elif pygem_prms.option_refreezing == 'Woodward': + 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 == 0: + self.surfacetype, self.firnline_idx = self._surfacetypebinsinitial(self.heights) + self.glac_bin_surfacetype_annual[:,year] = self.surfacetype + + # Off-glacier area and indices + if option_areaconstant == False: + self.offglac_bin_area_annual[:,year] = glacier_area_initial - glacier_area_t0 + offglac_idx = np.where(self.offglac_bin_area_annual[:,year] > 0)[0] + #TODO offglacier idx is not defined as off glacier, it's the idex with glacier area shrink + #print('****************** off glacier idx is*******************',offglac_idx) + + # Functions currently set up for monthly timestep + # only compute mass balance while glacier exists + if (pygem_prms.timestep == 'monthly'): + # if (pygem_prms.timestep == 'monthly') and (glac_idx_t0.shape[0] != 0): + + # AIR TEMPERATURE: Downscale the gcm temperature [deg C] to each bin + if pygem_prms.option_temp2bins == 1: + # 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[:,12*year:12*(year+1)] = (self.glacier_gcm_temp[12*year:12*(year+1)] + + self.glacier_gcm_lrgcm[12*year:12*(year+1)] * + (self.glacier_rgi_table.loc[pygem_prms.option_elev_ref_downscale] - self.glacier_gcm_elev) + + self.glacier_gcm_lrglac[12*year:12*(year+1)] * (heights - + self.glacier_rgi_table.loc[pygem_prms.option_elev_ref_downscale])[:, np.newaxis] + + self.modelprms['tbias']) + + # PRECIPITATION/ACCUMULATION: Downscale the precipitation (liquid and solid) to each bin + if pygem_prms.option_prec2bins == 1: + # Precipitation using precipitation factor and precipitation gradient + # P_bin = P_gcm * prec_factor * (1 + prec_grad * (z_bin - z_ref)) + bin_precsnow[:,12*year:12*(year+1)] = (self.glacier_gcm_prec[12*year:12*(year+1)] * + self.modelprms['kp'] * (1 + self.modelprms['precgrad'] * (heights - + self.glacier_rgi_table.loc[pygem_prms.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.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,12*year:12*(year+1)] = ( + bin_precsnow[glac_idx_75,12*year:12*(year+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 + 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.option_accumulation == 1: + # if temperature above threshold, then rain + (self.bin_prec[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold']]) = ( + bin_precsnow[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold']]) + # if temperature below threshold, then snow + (self.bin_acc[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] <= self.modelprms['tsnow_threshold']]) = ( + bin_precsnow[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] <= self.modelprms['tsnow_threshold']]) + elif pygem_prms.option_accumulation == 2: + # if temperature between min/max, then mix of snow/rain using linear relationship between min/max + self.bin_prec[:,12*year:12*(year+1)] = ( + (0.5 + (self.bin_temp[:,12*year:12*(year+1)] - + self.modelprms['tsnow_threshold']) / 2) * bin_precsnow[:,12*year:12*(year+1)]) + self.bin_acc[:,12*year:12*(year+1)] = ( + bin_precsnow[:,12*year:12*(year+1)] - self.bin_prec[:,12*year:12*(year+1)]) + # if temperature above maximum threshold, then all rain + (self.bin_prec[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold'] + 1]) = ( + bin_precsnow[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold'] + 1]) + (self.bin_acc[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] > self.modelprms['tsnow_threshold'] + 1]) = 0 + # if temperature below minimum threshold, then all snow + (self.bin_acc[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] <= self.modelprms['tsnow_threshold'] - 1]) = ( + bin_precsnow[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] <= self.modelprms['tsnow_threshold'] - 1]) + (self.bin_prec[:,12*year:12*(year+1)] + [self.bin_temp[:,12*year:12*(year+1)] <= 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 = 12*year + 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: - 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 - 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.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.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.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.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]]) - - # COMPUTE REFREEZING - TAP INTO "COLD RESERVOIR" or potential refreezing + 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.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.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]) + # 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.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.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.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.rf_dsc + + if pygem_prms.option_rf_limit_meltsnow == 1: + bin_meltlimit = self.bin_meltsnow.copy() else: + bin_meltlimit = self.bin_melt.copy() - if self.debug_refreeze and gidx == gidx_debug and step < 12: - print('\nMonth ' + str(self.dates_table.loc[step,'month']), 'Computing refreeze') + # Debug lowest bin + if self.debug_refreeze: + gidx_debug = np.where(heights == heights[glac_idx_t0].min())[0] - # Refreezing over firn surface - if (self.surfacetype[gidx] == 2) or (self.surfacetype[gidx] == 3): - nlayers = pygem_prms.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.pp) / pygem_prms.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.rf_layers - 1: - smax = pygem_prms.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: + # 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.rf_meltcrit: if self.debug_refreeze and gidx == gidx_debug and step < 12: - print('calculating potential refreeze from ' + str(nlayers) + ' layers') + print('\nMonth ' + str(self.dates_table.loc[step,'month']), + 'Computing heat conduction') - 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.rf_dz / pygem_prms.Lh_rf / pygem_prms.density_water) - 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)) + # 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.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.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.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('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] + print('tl_rf:', ["{:.2f}".format(x) for x in self.tl_rf[:,gidx,step]]) + + # COMPUTE REFREEZING - TAP INTO "COLD RESERVOIR" or potential refreezing 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 - - # Record refreeze - self.bin_refreeze[gidx,step] = self.refr[gidx] - - 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)) - - elif pygem_prms.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[:,12*year:12*(year+1)], - self.dates_table.iloc[12*year:12*(year+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('Year ' + str(year) + ' Month ' + str(self.dates_table.loc[step,'month']), - 'Refreeze potential:', np.round(refreeze_potential[glac_idx_t0[0]],3), - 'Snow depth:', np.round(self.bin_snowpack[glac_idx_t0[0],step],2), - 'Snow melt:', np.round(self.bin_meltsnow[glac_idx_t0[0],step],2), - 'Rain:', np.round(self.bin_prec[glac_idx_t0[0],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.tolerance, step] = 0 - # update refreeze potential - refreeze_potential -= self.bin_refreeze[:,step] - refreeze_potential[abs(refreeze_potential) < pygem_prms.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.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]) - - # 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] = self.glac_bin_massbalclim[:,12*year:12*(year+1)].sum(1) - # Update surface type for each bin - self.surfacetype, firnline_idx = self._surfacetypebinsannual(self.surfacetype, - self.glac_bin_massbalclim_annual, year) - # Record binned glacier area - self.glac_bin_area_annual[:,year] = glacier_area_t0 - # Store glacier-wide results - self._convert_glacwide_results(year, glacier_area_t0, heights, fls=fls, fl_id=fl_id, - option_areaconstant=option_areaconstant) - -## if debug: -# debug_startyr = 57 -# debug_endyr = 61 -# if year > debug_startyr and year < debug_endyr: -# print('\n', year, 'glac_bin_massbalclim:', self.glac_bin_massbalclim[:,12*year:12*(year+1)].sum(1)) -# print('ice thickness:', icethickness_t0) -# print('heights:', heights[glac_idx_t0]) -## print('surface type present:', self.glac_bin_surfacetype_annual[12:20,year]) -## print('surface type updated:', self.surfacetype[12:20]) - - # Mass balance for each bin [m ice per second] - seconds_in_year = self.dayspermonth[12*year:12*(year+1)].sum() * 24 * 3600 - mb = (self.glac_bin_massbalclim[:,12*year:12*(year+1)].sum(1) - * pygem_prms.density_water / pygem_prms.density_ice / seconds_in_year) - - if self.inversion_filter: - mb = np.minimum.accumulate(mb) - - # Fill in non-glaciated areas - needed for OGGM dynamics to remove small ice flux into next bin - mb_filled = mb.copy() - if len(glac_idx_t0) > 3: - mb_max = np.max(mb[glac_idx_t0]) - mb_min = np.min(mb[glac_idx_t0]) - height_max = np.max(heights[glac_idx_t0]) - height_min = np.min(heights[glac_idx_t0]) - mb_grad = (mb_min - mb_max) / (height_max - height_min) - 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]) - mb_filled[(mb_filled==0) & (heights < height_max)] = mb_min - -# if year > debug_startyr and year < debug_endyr: -# print('mb_min:', mb_min) -# -# if year > debug_startyr and year < debug_endyr: -# import matplotlib.pyplot as plt -# plt.plot(mb_filled, heights, '.') -# plt.ylabel('Elevation') -# plt.xlabel('Mass balance (mwea)') -# plt.show() -# -# print('mb_filled:', mb_filled) - - return mb_filled + if self.debug_refreeze and gidx == gidx_debug and step < 12: + print('\nMonth ' + str(self.dates_table.loc[step,'month']), 'Computing refreeze') + + # Refreezing over firn surface + if (self.surfacetype[gidx] == 2) or (self.surfacetype[gidx] == 3): + nlayers = pygem_prms.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.pp) / pygem_prms.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.rf_layers - 1: + smax = pygem_prms.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.rf_dz / pygem_prms.Lh_rf / pygem_prms.density_water) + 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 + + # Record refreeze + self.bin_refreeze[gidx,step] = self.refr[gidx] + + 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)) + + elif pygem_prms.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[:,12*year:12*(year+1)], + self.dates_table.iloc[12*year:12*(year+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('Year ' + str(year) + ' Month ' + str(self.dates_table.loc[step,'month']), + 'Refreeze potential:', np.round(refreeze_potential[glac_idx_t0[0]],3), + 'Snow depth:', np.round(self.bin_snowpack[glac_idx_t0[0],step],2), + 'Snow melt:', np.round(self.bin_meltsnow[glac_idx_t0[0],step],2), + 'Rain:', np.round(self.bin_prec[glac_idx_t0[0],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.tolerance, step] = 0 + # update refreeze potential + refreeze_potential -= self.bin_refreeze[:,step] + refreeze_potential[abs(refreeze_potential) < pygem_prms.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.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]) + + # 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] = self.glac_bin_massbalclim[:,12*year:12*(year+1)].sum(1) + # Update surface type for each bin + self.surfacetype, firnline_idx = self._surfacetypebinsannual(self.surfacetype, + self.glac_bin_massbalclim_annual, year) + # Record binned glacier area + self.glac_bin_area_annual[:,year] = glacier_area_t0 + # Store glacier-wide results + self._convert_glacwide_results(year, glacier_area_t0, heights, fls=fls, fl_id=fl_id, + option_areaconstant=option_areaconstant) + + ## if debug: + # debug_startyr = 57 + # debug_endyr = 61 + # if year > debug_startyr and year < debug_endyr: + # print('\n', year, 'glac_bin_massbalclim:', self.glac_bin_massbalclim[:,12*year:12*(year+1)].sum(1)) + # print('ice thickness:', icethickness_t0) + # print('heights:', heights[glac_idx_t0]) + ## print('surface type present:', self.glac_bin_surfacetype_annual[12:20,year]) + ## print('surface type updated:', self.surfacetype[12:20]) + + # Mass balance for each bin [m ice per second] + if year_month is None: + seconds_in_year = self.dayspermonth[12*year:12*(year+1)].sum() * 24 * 3600 + mb = (self.glac_bin_massbalclim[:,12*year:12*(year+1)].sum(1) + * pygem_prms.density_water / pygem_prms.density_ice / seconds_in_year) + else: + print("year in get_annul_mb is:",year) + print("year_month in get_annual_mb is:",year_month) + seconds_in_month = self.dayspermonth[12*year+int(year_month*12)]* 24 * 3600 + print("seconds_in_month",seconds_in_month) + mb = (self.glac_bin_massbalclim[:,12*year+int(year_month*12)] + * pygem_prms.density_water / pygem_prms.density_ice + /seconds_in_month) + print("index for mb") + print("Step (float)",12*year+int(year_month*12)) + print("Month",year_month) + print("Year (int)",year) + + if self.inversion_filter: + mb = np.minimum.accumulate(mb) + + # Fill in non-glaciated areas - needed for OGGM dynamics to remove small ice flux into next bin + mb_filled = mb.copy() + if len(glac_idx_t0) > 3: + mb_max = np.max(mb[glac_idx_t0]) + mb_min = np.min(mb[glac_idx_t0]) + height_max = np.max(heights[glac_idx_t0]) + height_min = np.min(heights[glac_idx_t0]) + mb_grad = (mb_min - mb_max) / (height_max - height_min) + 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]) + mb_filled[(mb_filled==0) & (heights < height_max)] = mb_min + + # if year > debug_startyr and year < debug_endyr: + # print('mb_min:', mb_min) + # + # if year > debug_startyr and year < debug_endyr: + # import matplotlib.pyplot as plt + # plt.plot(mb_filled, heights, '.') + # plt.ylabel('Elevation') + # plt.xlabel('Mass balance (mwea)') + # plt.show() + # + # print('mb_filled:', mb_filled) + print("**************** get_annual_mb end ****************") + return mb_filled + except: + print(traceback.format_exc()) #%% def _convert_glacwide_results(self, year, glacier_area, heights, @@ -734,6 +806,10 @@ def _convert_glacwide_results(self, year, glacier_area, heights, self.glac_wide_massbaltotal[12*year:12*(year+1)] = ( self.glac_wide_acc[12*year:12*(year+1)] + self.glac_wide_refreeze[12*year:12*(year+1)] - self.glac_wide_melt[12*year:12*(year+1)] - self.glac_wide_frontalablation[12*year:12*(year+1)]) + # Glacier-wide climatic mass balance (m3 w.e.) + self.glac_wide_massbalclim[12*year:12*(year+1)] = ( + self.glac_wide_acc[12*year:12*(year+1)] + self.glac_wide_refreeze[12*year:12*(year+1)] + - self.glac_wide_melt[12*year:12*(year+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: @@ -812,7 +888,7 @@ def _convert_glacwide_results(self, year, glacier_area, heights, ).sum(0)) - def ensure_mass_conservation(self, diag): + def ensure_mass_conservation(self, diag, Dynamic_step_Monthly = True): """ Ensure mass conservation that may result from using OGGM's glacier dynamics model. This will be resolved on an annual basis, and since the glacier dynamics are updated annually, the melt and runoff will be adjusted on a @@ -825,12 +901,27 @@ def ensure_mass_conservation(self, diag): Note: other dynamical models (e.g., mass redistribution curves, volume-length-area scaling) are based on the 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. + + Parameters + Dynamic_Step_Monthly : bool + if True, the dynamic step is monthly, the volume is monthly, need to be calculated as annual """ # Compute difference between volume change vol_change_annual_mbmod = (self.glac_wide_massbaltotal.reshape(-1,12).sum(1) * pygem_prms.density_water / pygem_prms.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] + #%% Dynamic running step + #TODO if the dynamic step is monthly, the volume is monthly, need to be calculated as annual + + if Dynamic_step_Monthly : + volume_m3_month_ice = diag.volume_m3.values[1:] - diag.volume_m3.values[0:-1] + volume_m3_annual_ice = volume_m3_month_ice.reshape(-1,12).sum(1) + vol_change_annual_diag = volume_m3_annual_ice + print("diag.volume_m3 :",diag.volume_m3) + print("vol_change_annual_diag :",vol_change_annual_diag) + else: + vol_change_annual_diag[0:diag.volume_m3.values[1:].shape[0]] = diag.volume_m3.values[1:] - diag.volume_m3.values[:-1] + #%% vol_change_annual_dif = vol_change_annual_diag - vol_change_annual_mbmod # Reduce glacier melt by the difference diff --git a/pygem/oggm_compat.py b/pygem/oggm_compat.py index 62e4a3d8..7cd94367 100755 --- a/pygem/oggm_compat.py +++ b/pygem/oggm_compat.py @@ -1,8 +1,10 @@ """ PYGEM-OGGGM COMPATIBILITY FUNCTIONS """ # External libraries +import os import numpy as np import pandas as pd import netCDF4 +import traceback # Local libraries import pygem_input as pygem_prms from oggm import cfg, utils @@ -12,6 +14,9 @@ from oggm.core.massbalance import MassBalanceModel #from oggm.shop import rgitopo from pygem.shop import debris, mbdata, icethickness +# add the millan22 thickness +from oggm.shop import millan22 + class CompatGlacDir: def __init__(self, rgiid): @@ -50,7 +55,7 @@ def single_flowline_glacier_directory(rgi_id, reset=pygem_prms.overwrite_gdirs, cfg.PARAMS['use_multiprocessing'] = False # Avoid erroneous glaciers (e.g., Centerlines too short or other issues) - cfg.PARAMS['continue_on_error'] = True + cfg.PARAMS['continue_on_error'] = False # Has internet cfg.PARAMS['has_internet'] = has_internet @@ -69,6 +74,8 @@ def single_flowline_glacier_directory(rgi_id, reset=pygem_prms.overwrite_gdirs, if not reset: try: gdir = utils.GlacierDirectory(rgi_id) + #gdir.read_pickle('inversion_flowlines', FILESUFFIX=) + #TODO Add suffix to the read_pickle function for diffferent parameters gdir.read_pickle('inversion_flowlines') # If the above works the directory is already processed, return return gdir @@ -79,36 +86,43 @@ def single_flowline_glacier_directory(rgi_id, reset=pygem_prms.overwrite_gdirs, process_gdir = True if process_gdir: - # Start after the prepro task level - base_url = pygem_prms.oggm_base_url - - cfg.PARAMS['has_internet'] = pygem_prms.has_internet - gdirs = workflow.init_glacier_directories([rgi_id], from_prepro_level=2, prepro_border=cfg.PARAMS['border'], - prepro_base_url=base_url, prepro_rgi_version='62') - - # Compute all the stuff - list_tasks = [ - # Consensus ice thickness - icethickness.consensus_gridded, - # Mass balance data - mbdata.mb_df_to_gdir] + try: + # Start after the prepro task level + base_url = pygem_prms.oggm_base_url + + cfg.PARAMS['has_internet'] = pygem_prms.has_internet + gdir = workflow.init_glacier_directories([rgi_id], from_prepro_level=2, prepro_border=cfg.PARAMS['border'], + prepro_base_url=base_url, prepro_rgi_version='62')[0] + # add the millan22 thickness and velocity + #workflow.execute_entity_task(millan22.velocity_to_gdir, gdirs) + #workflow.execute_entity_task(millan22.thickness_to_gdir, gdirs) + + # go through shop tasks to process auxiliary datasets to gdir if necessary + # consensus glacier mass + if not os.path.isfile(gdir.get_filepath('consensus_mass')): + workflow.execute_entity_task(icethickness.consensus_gridded, gdir) + # mass balance calibration data + if not os.path.isfile(gdir.get_filepath('mb_obs')): + workflow.execute_entity_task(mbdata.mb_df_to_gdir, gdir) - # Debris tasks - if pygem_prms.include_debris: - list_tasks.append(debris.debris_to_gdir) - list_tasks.append(debris.debris_binned) - - for task in list_tasks: - workflow.execute_entity_task(task, gdirs) - - gdir = gdirs[0] + # Debris tasks + # debris thickness and melt enhancement factors + if pygem_prms.include_debris: + if not os.path.isfile(gdir.get_filepath('debris_ed')) or os.path.isfile(gdir.get_filepath('debris_hd')): + workflow.execute_entity_task(debris.debris_to_gdir, gdir) + workflow.execute_entity_task(debris.debris_binned, gdir) + + gdir = gdir + except: + print("Somenthing wrong with the single_flowline_glacier_directory") + print(traceback.format_exc()) return gdir def single_flowline_glacier_directory_with_calving(rgi_id, reset=pygem_prms.overwrite_gdirs, - prepro_border=pygem_prms.oggm_border, k_calving=1, + prepro_border=pygem_prms.oggm_border, k_calving_str="", logging_level=pygem_prms.logging_level, has_internet=pygem_prms.has_internet): """Prepare a GlacierDirectory for PyGEM (single flowline to start with) @@ -143,7 +157,7 @@ def single_flowline_glacier_directory_with_calving(rgi_id, reset=pygem_prms.over cfg.PARAMS['use_multiprocessing'] = False # Avoid erroneous glaciers (e.g., Centerlines too short or other issues) - cfg.PARAMS['continue_on_error'] = True + cfg.PARAMS['continue_on_error'] = False # Has internet cfg.PARAMS['has_internet'] = has_internet @@ -156,7 +170,7 @@ def single_flowline_glacier_directory_with_calving(rgi_id, reset=pygem_prms.over cfg.PARAMS['dl_verify'] = True cfg.PARAMS['use_multiple_flowlines'] = False # temporary directory for testing (deleted on computer restart) - cfg.PATHS['working_dir'] = pygem_prms.oggm_gdir_fp + cfg.PATHS['working_dir'] = pygem_prms.oggm_gdir_fp+k_calving_str # Check if folder is already processed if not reset: @@ -175,24 +189,29 @@ def single_flowline_glacier_directory_with_calving(rgi_id, reset=pygem_prms.over # Start after the prepro task level base_url = pygem_prms.oggm_base_url - gdirs = workflow.init_glacier_directories([rgi_id], from_prepro_level=2, prepro_border=cfg.PARAMS['border'], - prepro_base_url=base_url, prepro_rgi_version='62') + gdir = workflow.init_glacier_directories([rgi_id], from_prepro_level=2, prepro_border=cfg.PARAMS['border'], + prepro_base_url=base_url, prepro_rgi_version='62')[0] + + #TODO: check if the glacier is tidewater, here we assume them to be tidewater whose calving calibarion is being done + gdir.is_tidewater = True - if not gdirs[0].is_tidewater: + if not gdir.is_tidewater: raise ValueError(f'{rgi_id} is not tidewater!') - # Compute all the stuff - list_tasks = [ - # Consensus ice thickness - icethickness.consensus_gridded, - # Mass balance data - mbdata.mb_df_to_gdir] + # add the millan22 thickness and velocity + #workflow.execute_entity_task(millan22.velocity_to_gdir, gdirs) + #workflow.execute_entity_task(millan22.thickness_to_gdir, gdirs) + + # go through shop tasks to process auxiliary datasets to gdir if necessary + # consensus glacier mass + if not os.path.isfile(gdir.get_filepath('consensus_mass')): + workflow.execute_entity_task(icethickness.consensus_gridded, gdir) - for task in list_tasks: - # The order matters! - workflow.execute_entity_task(task, gdirs) + # mass balance calibration data (note facorrected kwarg) + if not os.path.isfile(gdir.get_filepath('mb_obs')): + workflow.execute_entity_task(mbdata.mb_df_to_gdir, gdir) - return gdirs[0] + return gdir def create_empty_glacier_directory(rgi_id): diff --git a/pygem/output.py b/pygem/output.py new file mode 100644 index 00000000..c499957f --- /dev/null +++ b/pygem/output.py @@ -0,0 +1,660 @@ +#!/usr/bin/env python3 +""" +Created on Sept 19 2023 +Updated Mar 29 2024 + +@author: btobers mrweathers drounce + +PyGEM classes and subclasses for model output datasets + +For glacier simulations: +The two main parent classes are single_glacier(object) and compiled_regional(object) +Both of these have several subclasses which will inherit the necessary parent information + +""" +import pygem_input as pygem_prms +from dataclasses import dataclass +from scipy.stats import median_abs_deviation +from datetime import datetime +import numpy as np +import pandas as pd +import xarray as xr +import os, types, json, cftime, collections + +### single glacier output parent class ### +@dataclass +class single_glacier: + """ + Single glacier output dataset class for the Python Glacier Evolution Model. + """ + glacier_rgi_table : pd.DataFrame + dates_table : pd.DataFrame + pygem_version : float + gcm_name : str + scenario : str + realization : str + sim_iters : int + modelprms : dict + gcm_bc_startyear : int + gcm_startyear : int + gcm_endyear: int + + def __post_init__(self): + self.glac_values = np.array([self.glacier_rgi_table.name]) + self.glacier_str = '{0:0.5f}'.format(self.glacier_rgi_table['RGIId_float']) + self.reg_str = str(self.glacier_rgi_table.O1Region).zfill(2) + self.outdir = pygem_prms.output_sim_fp + self.set_fn() + self.set_time_vals() + self.model_params_record() + self.init_dicts() + + def set_fn(self, fn=None): + if fn: + self.outfn = fn + else: + self.outfn = self.glacier_str + '_' + self.gcm_name + '_' + if self.scenario: + self.outfn += f'{self.scenario}_' + if self.realization: + self.outfn += f'{self.realization}_' + if pygem_prms.option_calibration: + self.outfn += f'{pygem_prms.option_calibration}_' + else: + self.outfn += f'kp{self.modelprms["kp"]}_ddfsnow{self.modelprms["ddfsnow"]}_tbias{self.modelprms["tbias"]}_' + if self.gcm_name not in ['ERA-Interim', 'ERA5', 'COAWST']: + self.outfn += f'ba{pygem_prms.option_bias_adjustment}_' + yr0 = self.gcm_bc_startyear + else: + self.outfn += 'ba0_' + yr0 = self.gcm_startyear + if pygem_prms.option_calibration: + self.outfn += 'SETS_' + self.outfn += f'{yr0}_' + self.outfn += f'{self.gcm_endyear}_' + + def get_fn(self): + return self.outfn + + def set_time_vals(self): + if pygem_prms.gcm_wateryear == 'hydro': + self.year_type = 'water year' + self.annual_columns = np.unique(self.dates_table['wateryear'].values)[0:int(self.dates_table.shape[0]/12)] + elif pygem_prms.gcm_wateryear == 'calendar': + self.year_type = 'calendar year' + self.annual_columns = np.unique(self.dates_table['year'].values)[0:int(self.dates_table.shape[0]/12)] + elif pygem_prms.gcm_wateryear == 'custom': + self.year_type = 'custom year' + self.time_values = self.dates_table.loc[pygem_prms.gcm_spinupyears*12:self.dates_table.shape[0]+1,'date'].tolist() + self.time_values = [cftime.DatetimeNoLeap(x.year, x.month, x.day) for x in self.time_values] + # append additional year to self.year_values to account for mass and area at end of period + self.year_values = self.annual_columns[pygem_prms.gcm_spinupyears:self.annual_columns.shape[0]] + self.year_values = np.concatenate((self.year_values, np.array([self.annual_columns[-1] + 1]))) + + def model_params_record(self): + substrings = ['user_info', 'rgi', 'glac_n', 'fp', 'fn', 'filepath', 'directory','url','logging','overwrite','hugonnet'] # substrings to look for in pygem_prms that don't necessarily need to store + # get all locally defined variables from the pygem_prms, excluding imports, functions, and classes + self.mdl_params_dict = { + var: value + for var, value in vars(pygem_prms).items() + if not var.startswith('__') and + not isinstance(value, (types.ModuleType, types.FunctionType, type)) and + not any(substring.lower() in var.lower() for substring in substrings) + } + # overwrite variables that are possibly different from pygem_input + self.mdl_params_dict['gcm_bc_startyear'] = self.gcm_bc_startyear + self.mdl_params_dict['gcm_startyear'] = self.gcm_startyear + self.mdl_params_dict['gcm_endyear'] = self.gcm_endyear + self.mdl_params_dict['gcm_name'] = self.gcm_name + self.mdl_params_dict['realization'] = self.realization + self.mdl_params_dict['scenario'] = self.scenario + + def init_dicts(self): + # Boilerplate coordinate and attribute dictionaries - these will be the same for both glacier-wide and binned outputs + self.output_coords_dict = collections.OrderedDict() + self.output_coords_dict['RGIId'] = collections.OrderedDict([('glac', self.glac_values)]) + self.output_coords_dict['CenLon'] = collections.OrderedDict([('glac', self.glac_values)]) + self.output_coords_dict['CenLat'] = collections.OrderedDict([('glac', self.glac_values)]) + self.output_coords_dict['O1Region'] = collections.OrderedDict([('glac', self.glac_values)]) + self.output_coords_dict['O2Region'] = collections.OrderedDict([('glac', self.glac_values)]) + self.output_coords_dict['Area'] = collections.OrderedDict([('glac', self.glac_values)]) + self.output_attrs_dict = { + 'time': { + 'long_name': 'time', + 'year_type':self.year_type, + 'comment':'start of the month'}, + 'glac': { + 'long_name': 'glacier index', + 'comment': 'glacier index referring to glaciers properties and model results'}, + 'year': { + 'long_name': 'years', + 'year_type': self.year_type, + 'comment': 'years referring to the start of each year'}, + 'RGIId': { + 'long_name': 'Randolph Glacier Inventory ID', + 'comment': 'RGIv6.0'}, + 'CenLon': { + 'long_name': 'center longitude', + 'units': 'degrees E', + 'comment': 'value from RGIv6.0'}, + 'CenLat': { + 'long_name': 'center latitude', + 'units': 'degrees N', + 'comment': 'value from RGIv6.0'}, + 'O1Region': { + 'long_name': 'RGI order 1 region', + 'comment': 'value from RGIv6.0'}, + 'O2Region': { + 'long_name': 'RGI order 2 region', + 'comment': 'value from RGIv6.0'}, + 'Area': { + 'long_name': 'glacier area', + 'units': 'm2', + 'comment': 'value from RGIv6.0'} + } + def create_xr_ds(self): + # Add variables to empty dataset and merge together + count_vn = 0 + self.encoding = {} + for vn in self.output_coords_dict.keys(): + count_vn += 1 + empty_holder = np.zeros([len(self.output_coords_dict[vn][i]) for i in list(self.output_coords_dict[vn].keys())]) + output_xr_ds_ = xr.Dataset({vn: (list(self.output_coords_dict[vn].keys()), empty_holder)}, + coords=self.output_coords_dict[vn]) + # Merge datasets of stats into one output + if count_vn == 1: + self.output_xr_ds = output_xr_ds_ + else: + self.output_xr_ds = xr.merge((self.output_xr_ds, output_xr_ds_)) + noencoding_vn = ['RGIId'] + # Add attributes + for vn in self.output_xr_ds.variables: + try: + self.output_xr_ds[vn].attrs = self.output_attrs_dict[vn] + except: + pass + # Encoding (specify _FillValue, offsets, etc.) + + if vn not in noencoding_vn: + self.encoding[vn] = {'_FillValue': None, + 'zlib':True, + 'complevel':9 + } + self.output_xr_ds['RGIId'].values = np.array([self.glacier_rgi_table.loc['RGIId']]) + self.output_xr_ds['CenLon'].values = np.array([self.glacier_rgi_table.CenLon]) + self.output_xr_ds['CenLat'].values = np.array([self.glacier_rgi_table.CenLat]) + self.output_xr_ds['O1Region'].values = np.array([self.glacier_rgi_table.O1Region]) + self.output_xr_ds['O2Region'].values = np.array([self.glacier_rgi_table.O2Region]) + self.output_xr_ds['Area'].values = np.array([self.glacier_rgi_table.Area * 1e6]) + + self.output_xr_ds.attrs = {'source': f'PyGEMv{self.pygem_version}', + 'institution': pygem_prms.user_info['institution'], + 'history': f'Created by {pygem_prms.user_info["name"]} ({pygem_prms.user_info["email"]}) on ' + datetime.today().strftime('%Y-%m-%d'), + 'references': 'doi:10.3389/feart.2019.00331 and doi:10.1017/jog.2019.91', + 'model_parameters':json.dumps(self.mdl_params_dict)} + + def get_xr_ds(self): + return self.output_xr_ds + + def save_xr_ds(self, netcdf_fn): + # export netcdf + self.output_xr_ds.to_netcdf(self.outdir + netcdf_fn, encoding=self.encoding) + # close datasets + self.output_xr_ds.close() + + +@dataclass +class glacierwide_stats(single_glacier): + """ + Single glacier-wide statistics dataset + """ + + def __post_init__(self): + super().__post_init__() # call parent class __post_init__ (get glacier values, time stamps, and instantiate output dictionaries that will form netcdf file output) + self.set_outdir() + self.update_dicts() # add required fields to output dictionary + + def set_outdir(self): + # prepare for export + self.outdir += self.reg_str + '/' + self.gcm_name + '/' + if self.gcm_name not in ['ERA-Interim', 'ERA5', 'COAWST']: + self.outdir += self.scenario + '/' + self.outdir += 'stats/' + # Create filepath if it does not exist + os.makedirs(self.outdir, exist_ok=True) + + def update_dicts(self): + # update coordinate and attribute dictionaries + self.output_coords_dict['glac_runoff_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_runoff_monthly'] = { + 'long_name': 'glacier-wide runoff', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'runoff from the glacier terminus, which moves over time'} + self.output_coords_dict['glac_area_annual'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_area_annual'] = { + 'long_name': 'glacier area', + 'units': 'm2', + 'temporal_resolution': 'annual', + 'comment': 'area at start of the year'} + self.output_coords_dict['glac_mass_annual'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_mass_annual'] = { + 'long_name': 'glacier mass', + 'units': 'kg', + 'temporal_resolution': 'annual', + 'comment': 'mass of ice based on area and ice thickness at start of the year'} + self.output_coords_dict['glac_mass_bsl_annual'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_mass_bsl_annual'] = { + 'long_name': 'glacier mass below sea level', + 'units': 'kg', + 'temporal_resolution': 'annual', + 'comment': 'mass of ice below sea level based on area and ice thickness at start of the year'} + self.output_coords_dict['glac_ELA_annual'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_ELA_annual'] = { + 'long_name': 'annual equilibrium line altitude above mean sea level', + 'units': 'm', + '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([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_runoff_monthly'] = { + 'long_name': 'off-glacier-wide runoff', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'off-glacier runoff from area where glacier no longer exists'} + if self.sim_iters > 1: + self.output_coords_dict['glac_runoff_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_runoff_monthly_mad'] = { + 'long_name': 'glacier-wide runoff median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'runoff from the glacier terminus, which moves over time'} + self.output_coords_dict['glac_area_annual_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_area_annual_mad'] = { + 'long_name': 'glacier area median absolute deviation', + 'units': 'm2', + 'temporal_resolution': 'annual', + 'comment': 'area at start of the year'} + self.output_coords_dict['glac_mass_annual_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_mass_annual_mad'] = { + 'long_name': 'glacier mass median absolute deviation', + 'units': 'kg', + 'temporal_resolution': 'annual', + 'comment': 'mass of ice based on area and ice thickness at start of the year'} + self.output_coords_dict['glac_mass_bsl_annual_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_mass_bsl_annual_mad'] = { + 'long_name': 'glacier mass below sea level median absolute deviation', + 'units': 'kg', + 'temporal_resolution': 'annual', + 'comment': 'mass of ice below sea level based on area and ice thickness at start of the year'} + self.output_coords_dict['glac_ELA_annual_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_ELA_annual_mad'] = { + 'long_name': 'annual equilibrium line altitude above mean sea level median absolute deviation', + 'units': 'm', + '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([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_runoff_monthly_mad'] = { + 'long_name': 'off-glacier-wide runoff median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'off-glacier runoff from area where glacier no longer exists'} + + if pygem_prms.export_extra_vars: + self.output_coords_dict['glac_prec_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_prec_monthly'] = { + 'long_name': 'glacier-wide precipitation (liquid)', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'only the liquid precipitation, solid precipitation excluded'} + self.output_coords_dict['glac_temp_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_temp_monthly'] = { + 'standard_name': 'air_temperature', + 'long_name': 'glacier-wide mean air temperature', + 'units': 'K', + 'temporal_resolution': 'monthly', + '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([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_acc_monthly'] = { + 'long_name': 'glacier-wide accumulation, in water equivalent', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'only the solid precipitation'} + self.output_coords_dict['glac_refreeze_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_refreeze_monthly'] = { + 'long_name': 'glacier-wide refreeze, in water equivalent', + 'units': 'm3', + 'temporal_resolution': 'monthly'} + self.output_coords_dict['glac_melt_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_melt_monthly'] = { + 'long_name': 'glacier-wide melt, in water equivalent', + 'units': 'm3', + 'temporal_resolution': 'monthly'} + self.output_coords_dict['glac_frontalablation_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_frontalablation_monthly'] = { + 'long_name': 'glacier-wide frontal ablation, in water equivalent', + 'units': 'm3', + 'temporal_resolution': 'monthly', + '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([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_massbaltotal_monthly'] = { + 'long_name': 'glacier-wide total mass balance, in water equivalent', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'total mass balance is the sum of the climatic mass balance and frontal ablation'} + self.output_coords_dict['glac_snowline_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_snowline_monthly'] = { + 'long_name': 'transient snowline altitude above mean sea level', + 'units': 'm', + 'temporal_resolution': 'monthly', + 'comment': 'transient snowline is altitude separating snow from ice/firn'} + self.output_coords_dict['glac_mass_change_ignored_annual'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_mass_change_ignored_annual'] = { + 'long_name': 'glacier mass change ignored', + 'units': 'kg', + 'temporal_resolution': 'annual', + 'comment': 'glacier mass change ignored due to flux divergence'} + self.output_coords_dict['offglac_prec_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_prec_monthly'] = { + 'long_name': 'off-glacier-wide precipitation (liquid)', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'only the liquid precipitation, solid precipitation excluded'} + self.output_coords_dict['offglac_refreeze_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_refreeze_monthly'] = { + 'long_name': 'off-glacier-wide refreeze, in water equivalent', + 'units': 'm3', + 'temporal_resolution': 'monthly'} + self.output_coords_dict['offglac_melt_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_melt_monthly'] = { + 'long_name': 'off-glacier-wide melt, in water equivalent', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'only melt of snow and refreeze since off-glacier'} + self.output_coords_dict['offglac_snowpack_monthly'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_snowpack_monthly'] = { + 'long_name': 'off-glacier-wide snowpack, in water equivalent', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'snow remaining accounting for new accumulation, melt, and refreeze'} + + if self.sim_iters > 1: + self.output_coords_dict['glac_prec_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_prec_monthly_mad'] = { + 'long_name': 'glacier-wide precipitation (liquid) median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'only the liquid precipitation, solid precipitation excluded'} + self.output_coords_dict['glac_temp_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_temp_monthly_mad'] = { + 'standard_name': 'air_temperature', + 'long_name': 'glacier-wide mean air temperature median absolute deviation', + 'units': 'K', + 'temporal_resolution': 'monthly', + '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([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_acc_monthly_mad'] = { + 'long_name': 'glacier-wide accumulation, in water equivalent, median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'only the solid precipitation'} + self.output_coords_dict['glac_refreeze_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_refreeze_monthly_mad'] = { + 'long_name': 'glacier-wide refreeze, in water equivalent, median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly'} + self.output_coords_dict['glac_melt_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_melt_monthly_mad'] = { + 'long_name': 'glacier-wide melt, in water equivalent, median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly'} + self.output_coords_dict['glac_frontalablation_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_frontalablation_monthly_mad'] = { + 'long_name': 'glacier-wide frontal ablation, in water equivalent, median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly', + '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([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_massbaltotal_monthly_mad'] = { + 'long_name': 'glacier-wide total mass balance, in water equivalent, median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly', + '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([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['glac_snowline_monthly_mad'] = { + 'long_name': 'transient snowline above mean sea level median absolute deviation', + 'units': 'm', + 'temporal_resolution': 'monthly', + 'comment': 'transient snowline is altitude separating snow from ice/firn'} + self.output_coords_dict['glac_mass_change_ignored_annual_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('year', self.year_values)]) + self.output_attrs_dict['glac_mass_change_ignored_annual_mad'] = { + 'long_name': 'glacier mass change ignored median absolute deviation', + 'units': 'kg', + 'temporal_resolution': 'annual', + 'comment': 'glacier mass change ignored due to flux divergence'} + self.output_coords_dict['offglac_prec_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_prec_monthly_mad'] = { + 'long_name': 'off-glacier-wide precipitation (liquid) median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'only the liquid precipitation, solid precipitation excluded'} + self.output_coords_dict['offglac_refreeze_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_refreeze_monthly_mad'] = { + 'long_name': 'off-glacier-wide refreeze, in water equivalent, median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly'} + self.output_coords_dict['offglac_melt_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_melt_monthly_mad'] = { + 'long_name': 'off-glacier-wide melt, in water equivalent, median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'only melt of snow and refreeze since off-glacier'} + self.output_coords_dict['offglac_snowpack_monthly_mad'] = collections.OrderedDict([('glac', self.glac_values), + ('time', self.time_values)]) + self.output_attrs_dict['offglac_snowpack_monthly_mad'] = { + 'long_name': 'off-glacier-wide snowpack, in water equivalent, median absolute deviation', + 'units': 'm3', + 'temporal_resolution': 'monthly', + 'comment': 'snow remaining accounting for new accumulation, melt, and refreeze'} + + +@dataclass +class binned_stats(single_glacier): + """ + Single glacier binned dataset + """ + nbins : int + + def __post_init__(self): + super().__post_init__() # call parent class __post_init__ (get glacier values, time stamps, and instantiate output dictionaries that will form netcdf file output) + self.bin_values = np.arange(self.nbins) # bin indices + self.set_outdir() + self.update_dicts() # add required fields to output dictionary + + def set_outdir(self): + # prepare for export + self.outdir += self.reg_str + '/' + self.gcm_name + '/' + if self.gcm_name not in ['ERA-Interim', 'ERA5', 'COAWST']: + self.outdir += self.scenario + '/' + self.outdir += 'binned/' + # Create filepath if it does not exist + os.makedirs(self.outdir, exist_ok=True) + + def update_dicts(self): + # update coordinate and attribute dictionaries + self.output_coords_dict['bin_distance'] = collections.OrderedDict([('glac', self.glac_values), ('bin', self.bin_values)]) + self.output_attrs_dict['bin_distance'] = { + 'long_name': 'distance downglacier', + 'units': 'm', + 'comment': 'horizontal distance calculated from top of glacier moving downglacier'} + self.output_coords_dict['bin_surface_h_initial'] = collections.OrderedDict([('glac', self.glac_values), ('bin', self.bin_values)]) + self.output_attrs_dict['bin_surface_h_initial'] = { + 'long_name': 'initial binned surface elevation', + 'units': 'm above sea level'} + self.output_coords_dict['bin_mass_annual'] = ( + collections.OrderedDict([('glac', self.glac_values), ('bin', self.bin_values), ('year', self.year_values)])) + self.output_attrs_dict['bin_mass_annual'] = { + 'long_name': 'binned ice mass', + 'units': 'kg', + 'temporal_resolution': 'annual', + 'comment': 'binned ice mass at start of the year'} + self.output_coords_dict['bin_thick_annual'] = ( + collections.OrderedDict([('glac', self.glac_values), ('bin', self.bin_values), ('year', self.year_values)])) + self.output_attrs_dict['bin_thick_annual'] = { + 'long_name': 'binned ice thickness', + 'units': 'm', + 'temporal_resolution': 'annual', + 'comment': 'binned ice thickness at start of the year'} + self.output_coords_dict['bin_massbalclim_annual'] = ( + collections.OrderedDict([('glac', self.glac_values), ('bin', self.bin_values), ('year', self.year_values)])) + self.output_attrs_dict['bin_massbalclim_annual'] = { + 'long_name': 'binned climatic mass balance, in water equivalent', + 'units': 'm', + 'temporal_resolution': 'annual', + 'comment': 'climatic mass balance is computed before dynamics so can theoretically exceed ice thickness'}, + self.output_coords_dict['bin_massbalclim_monthly'] = ( + 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', + 'units': 'm', + 'temporal_resolution': 'monthly', + 'comment': 'monthly climatic mass balance from the PyGEM mass balance module'} + + if self.sim_iters > 1: + self.output_coords_dict['bin_mass_annual_mad'] = ( + collections.OrderedDict([('glac', self.glac_values), ('bin', self.bin_values), ('year', self.year_values)])) + self.output_attrs_dict['bin_mass_annual_mad'] = { + 'long_name': 'binned ice mass median absolute deviation', + 'units': 'kg', + 'temporal_resolution': 'annual', + 'comment': 'mass of ice based on area and ice thickness at start of the year'} + self.output_coords_dict['bin_thick_annual_mad'] = ( + collections.OrderedDict([('glac', self.glac_values), ('bin', self.bin_values), ('year', self.year_values)])) + self.output_attrs_dict['bin_thick_annual_mad'] = { + 'long_name': 'binned ice thickness median absolute deviation', + 'units': 'm', + 'temporal_resolution': 'annual', + 'comment': 'thickness of ice at start of the year'} + self.output_coords_dict['bin_massbalclim_annual_mad'] = ( + collections.OrderedDict([('glac', self.glac_values), ('bin', self.bin_values), ('year', self.year_values)])) + self.output_attrs_dict['bin_massbalclim_annual_mad'] = { + 'long_name': 'binned climatic mass balance, in water equivalent, median absolute deviation', + 'units': 'm', + 'temporal_resolution': 'annual', + 'comment': 'climatic mass balance is computed before dynamics so can theoretically exceed ice thickness'} + + + +### compiled regional output parent class ### +@dataclass +class compiled_regional: + """ + Compiled regional output dataset for the Python Glacier Evolution Model. + """ + +@dataclass +class regional_annual_mass(compiled_regional): + """ + compiled regional annual mass + """ + +@dataclass +class regional_annual_area(compiled_regional): + """ + compiled regional annual area + """ + +@dataclass +class regional_monthly_runoff(compiled_regional): + """ + compiled regional monthly runoff + """ + +@dataclass +class regional_monthly_massbal(compiled_regional): + """ + compiled regional monthly climatic mass balance + """ + + +def calc_stats_array(data, stats_cns=['median', 'mad']): + """ + Calculate stats for a given variable + + Parameters + ---------- + vn : str + variable name + ds : xarray dataset + dataset of output with all ensemble simulations + + Returns + ------- + stats : np.array + Statistics related to a given variable + """ + stats = None + if 'mean' in stats_cns: + if stats is None: + stats = np.nanmean(data,axis=1)[:,np.newaxis] + if 'std' in stats_cns: + stats = np.append(stats, np.nanstd(data,axis=1)[:,np.newaxis], axis=1) + if '2.5%' in stats_cns: + stats = np.append(stats, np.nanpercentile(data, 2.5, axis=1)[:,np.newaxis], axis=1) + if '25%' in stats_cns: + stats = np.append(stats, np.nanpercentile(data, 25, axis=1)[:,np.newaxis], axis=1) + if 'median' in stats_cns: + if stats is None: + stats = np.nanmedian(data, axis=1)[:,np.newaxis] + else: + stats = np.append(stats, np.nanmedian(data, axis=1)[:,np.newaxis], axis=1) + if '75%' in stats_cns: + stats = np.append(stats, np.nanpercentile(data, 75, axis=1)[:,np.newaxis], axis=1) + if '97.5%' in stats_cns: + stats = np.append(stats, np.nanpercentile(data, 97.5, axis=1)[:,np.newaxis], axis=1) + if 'mad' in stats_cns: + stats = np.append(stats, median_abs_deviation(data, axis=1, nan_policy='omit')[:,np.newaxis], axis=1) + return stats \ No newline at end of file diff --git a/pygem/pygem_modelsetup.py b/pygem/pygem_modelsetup.py index 7802976a..089a6b50 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -6,6 +6,7 @@ import pandas as pd import numpy as np from datetime import datetime +import traceback # Local libraries import pygem_input as pygem_prms @@ -56,7 +57,7 @@ def datesmodelrun(startyear=pygem_prms.ref_startyear, endyear=pygem_prms.ref_end if pygem_prms.timestep == 'monthly': # Automatically generate dates from start date to end data using a monthly frequency (MS), which generates # monthly data using the 1st of each month' - dates_table = pd.DataFrame({'date' : pd.date_range(startdate, enddate, freq='MS', unit='s')}) + dates_table = pd.DataFrame({'date' : pd.date_range(startdate, enddate, freq='MS')}) # 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 @@ -264,7 +265,8 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' indexname=pygem_prms.indexname, include_landterm=True,include_laketerm=True,include_tidewater=True, glac_no_skip=pygem_prms.glac_no_skip, - min_glac_area_km2=0): + min_glac_area_km2=0, + debug=False): """ Select all glaciers to be used in the model run according to the regions and glacier numbers defined by the RGI glacier inventory. This function returns the rgi table associated with all of these glaciers. @@ -297,32 +299,42 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' # Create an empty dataframe rgi_regionsO1 = sorted(rgi_regionsO1) + print('glc_no',glac_no) + #print('rfp',os.listdir(rgi_fp)) + print('rgi_glc_num',rgi_glac_number) glacier_table = pd.DataFrame() for region in rgi_regionsO1: - + print(region) if glac_no is not None: rgi_glac_number = glac_no_byregion[region] - + print("rgi_glac_number:",rgi_glac_number) # if len(rgi_glac_number) < 50: for i in os.listdir(rgi_fp): if i.startswith(str(region).zfill(2)) and i.endswith('.csv'): rgi_fn = i + #print('regs:',i) try: csv_regionO1 = pd.read_csv(rgi_fp + rgi_fn) except: csv_regionO1 = pd.read_csv(rgi_fp + rgi_fn, encoding='latin1') + print('..latin1 encoding..') + #print(traceback.format_exc()) # Populate glacer_table with the glaciers of interest if rgi_regionsO2 == 'all' and rgi_glac_number == 'all': - print("All glaciers within region(s) %s are included in this model run." % (region)) + if debug: + print("All glaciers within region(s) %s are included in this model run." % (region)) if glacier_table.empty: glacier_table = csv_regionO1 + print('..1...',glacier_table.shape[0]) else: glacier_table = pd.concat([glacier_table, csv_regionO1], axis=0) + print('..2...',glacier_table.shape[0]) elif rgi_regionsO2 != 'all' and rgi_glac_number == 'all': - print("All glaciers within subregion(s) %s in region %s are included in this model run." % - (rgi_regionsO2, region)) + if debug: + print("All glaciers within subregion(s) %s in region %s are included in this model run." % + (rgi_regionsO2, region)) for regionO2 in rgi_regionsO2: if glacier_table.empty: glacier_table = csv_regionO1.loc[csv_regionO1['O2Region'] == regionO2] @@ -345,8 +357,9 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' else: glacier_table = (pd.concat([glacier_table, csv_regionO1.loc[rgi_idx]], axis=0)) - + print("glacier_table:",glacier_table) glacier_table = glacier_table.copy() + print('num glc in the region:',glacier_table.shape[0]) # reset the index so that it is in sequential order (0, 1, 2, etc.) glacier_table.reset_index(inplace=True) # drop connectivity 2 for Greenland and Antarctica @@ -393,6 +406,7 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' termtype_values.append(5) if include_laketerm: termtype_values.append(2) + #print(glacier.columns.tolist()) glacier_table = glacier_table.loc[glacier_table['TermType'].isin(termtype_values)] glacier_table.reset_index(inplace=True, drop=True) # Glacier number with no trailing zeros @@ -494,7 +508,7 @@ def split_list(lst, n=1, option_ordered=1, group_thousands=False): merged = [item for sublist in lst_batches for item in sublist if item[:5]==s] lst_batches_th.append(merged) # ensure that number of batches doesn't exceed original number - while len(lst_batches_th) > len(lst_batches): + while len(lst_batches_th) > n: # move shortest batch to next shortest batch lengths = np.asarray([len(batch) for batch in lst_batches_th]) sorted = lengths.argsort() diff --git a/pygem/shop/debris.py b/pygem/shop/debris.py index 817a8510..70299d39 100755 --- a/pygem/shop/debris.py +++ b/pygem/shop/debris.py @@ -114,7 +114,7 @@ def debris_to_gdir(gdir, debris_dir=pygem_prms.debris_fp, add_to_gridded=True, h @entity_task(log, writes=['inversion_flowlines']) -def debris_binned(gdir, ignore_debris=False, fl_str='inversion_flowlines'): +def debris_binned(gdir, ignore_debris=False, fl_str='inversion_flowlines',filesuffix=''): """Bin debris thickness and enhancement factors. Updates the 'inversion_flowlines' save file. @@ -123,10 +123,16 @@ def debris_binned(gdir, ignore_debris=False, fl_str='inversion_flowlines'): ---------- gdir : :py:class:`oggm.GlacierDirectory` where to write the data + ignore_debris : bool + option to ignore debris data + fl_str : str + filename of inversion flowlines + filesuffix : str + filesuffix for inversion flow """ # Nominal glaciers will throw error, so make sure inversion_flowlines exist try: - flowlines = gdir.read_pickle(fl_str) + flowlines = gdir.read_pickle(fl_str, filesuffix=filesuffix) fl = flowlines[0] assert len(flowlines) == 1, 'Error: binning debris only works for single flowlines at present' @@ -186,5 +192,5 @@ def debris_binned(gdir, ignore_debris=False, fl_str='inversion_flowlines'): fl.debris_ed = np.ones(nbins) # Overwrite pickle - gdir.write_pickle(flowlines, fl_str) + gdir.write_pickle(flowlines, fl_str, filesuffix=filesuffix) \ No newline at end of file diff --git a/pygem/shop/mbdata.py b/pygem/shop/mbdata.py index 39e3093d..e5e82aef 100755 --- a/pygem/shop/mbdata.py +++ b/pygem/shop/mbdata.py @@ -79,6 +79,7 @@ def mb_df_to_gdir(gdir, mb_dataset='Hugonnet2020'): mb_mwea_err = mb_df.loc[rgiid_idx, mberr_cn] assert mb_clim_cn in mb_df.columns, mb_clim_cn + ' not a column in mb_df' + #print("in mb_df_to_gdir, mb_clim_cn is :",mb_clim_cn) mb_clim_mwea = mb_df.loc[rgiid_idx, mb_clim_cn] mb_clim_mwea_err = mb_df.loc[rgiid_idx, mberr_clim_cn] @@ -111,6 +112,7 @@ def mb_df_to_gdir(gdir, mb_dataset='Hugonnet2020'): 'nyears': nyears} pkl_fn = gdir.get_filepath('mb_obs') + print("the path of 'mb_obs'in gdir is :",pkl_fn) with open(pkl_fn, 'wb') as f: pickle.dump(mbdata, f) diff --git a/pyproject.toml b/pyproject.toml index 1cd5a853..3d35c1e3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,6 +22,7 @@ classifiers = [ "License :: OSI Approved :: MIT License", "Operating System :: OS Independent", ] +dynamic = ["license", "dependencies"] [project.urls] "Homepage" = "https://github.com/drounce/PyGEM" \ No newline at end of file diff --git a/setup.py b/setup.py old mode 100755 new mode 100644 index 46cfcbad..0430d327 --- a/setup.py +++ b/setup.py @@ -1,19 +1,6 @@ -#!/usr/bin/env python3 -import setuptools -from setuptools import find_packages # or find_namespace_packages +#!/usr/bin/env python +from setuptools import setup if __name__ == "__main__": - setuptools.setup( - # ..., - # package_dir={'':'pygem'} - # ... - # packages= [] - # find_packages( - # All keyword arguments below are optional: - - # where='pygem', # '.' by default - # include=['mypackage*'], # ['*'] by default - # exclude=['mypackage.tests'], # empty by default - # ), - # ... - ) \ No newline at end of file + + setup() \ No newline at end of file