From 73abeac2c96a80698e124110f5d2e581a368574b Mon Sep 17 00:00:00 2001 From: ruitang yang Date: Wed, 15 Nov 2023 14:05:26 +0100 Subject: [PATCH 01/45] make date_range calls compatable with pandas <2 --- pygem/pygem_modelsetup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygem/pygem_modelsetup.py b/pygem/pygem_modelsetup.py index 7802976a..c09339b9 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -56,7 +56,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 From dc8cab15e9e87af472d89e12c1bb016384e2c30b Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 4 Dec 2023 12:17:06 +0100 Subject: [PATCH 02/45] add the sermeq calving function fa_sermeq_speed_law --- pygem/glacierdynamics.py | 295 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 295 insertions(+) diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index fc6d7324..d5127cb5 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -455,8 +455,297 @@ 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() + + #%% ----- Calving Law ---------- +def fa_sermeq_speed_law(model,last_above_wl, v_scaling=1, verbose=False, + tau0=150e3, 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=150e3, 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. + """ + 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 = tau0 + mu * N + else: # assume constant if not set + ty = tau0 + 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)) ** 1)) + # 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 + + #%% ----- FRONTAL ABLATION ----- def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, calving_k=None, debug=False ): @@ -521,6 +810,12 @@ def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, 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] + # Sermeq calving law + + + + + # 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] From ab32306bfdd514d320116cb49cb25f2ad2997251 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 4 Dec 2023 12:31:35 +0100 Subject: [PATCH 03/45] replace the claving law in _get_annual_frontalablation by sermeq law --- pygem/glacierdynamics.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index d5127cb5..cc4d5d49 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -455,8 +455,8 @@ 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() - - + + #%% ----- Calving Law ---------- def fa_sermeq_speed_law(model,last_above_wl, v_scaling=1, verbose=False, tau0=150e3, yield_type='constant', mu=0.01, @@ -809,13 +809,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,v_scaling = 1, verbose = False,tau0 = 150e3, + 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] From 7d67e72beca20134d12fa682a67655c0e0553c0b Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 4 Dec 2023 15:28:22 +0100 Subject: [PATCH 04/45] set the value of tau0 as the value of calving_k --- pygem/glacierdynamics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index cc4d5d49..e8ef0a96 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -812,7 +812,7 @@ def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, #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,v_scaling = 1, verbose = False,tau0 = 150e3, + s_fa = fa_sermeq_speed_law(self,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********************') From 08d088c2c99fb02e19b7b7e7eae9ba2ceabca71f Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Thu, 7 Dec 2023 09:30:31 +0100 Subject: [PATCH 05/45] scale the calving-k(yield strength) --- pygem/glacierdynamics.py | 1351 +++++++++++++++++++------------------ pygem/pygem_modelsetup.py | 10 +- 2 files changed, 690 insertions(+), 671 deletions(-) diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index e8ef0a96..6d583743 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -21,729 +21,740 @@ 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 -#%% -class MassRedistributionCurveModel(FlowlineModel): - """Glacier geometry updated using mass redistribution curves; also known as the "delta-h method" + @author: Ruitang Yang & Lizz Ultee - This uses mass redistribution curves from Huss et al. (2010) to update the glacier geometry - """ + Authors: Ruitang Yang & Lizz Ultee + Parameters + ---------- - def __init__(self, flowlines, mb_model=None, y0=0., - glen_a=None, fs=0., is_tidewater=False, water_level=0, -# calving_k=0, - inplace=False, - debug=True, - option_areaconstant=False, spinupyears=0, - constantarea_years=0, - **kwargs): - """ Instanciate the model. - - Parameters - ---------- - flowlines : list - the glacier flowlines - mb_model : MassBalanceModel - the mass-balance model - y0 : int - initial year of the simulation - inplace : bool - whether or not to make a copy of the flowline objects for the run - setting to True implies that your objects will be modified at run - time by the model (can help to spare memory) - is_tidewater: bool, default: False - use the very basic parameterization for tidewater glaciers - mb_elev_feedback : str, default: 'annual' - 'never', 'always', 'annual', or 'monthly': how often the - mass-balance should be recomputed from the mass balance model. - 'Never' is equivalent to 'annual' but without elevation feedback - at all (the heights are taken from the first call). - check_for_boundaries: bool, default: True - raise an error when the glacier grows bigger than the domain - boundaries - """ - super(MassRedistributionCurveModel, self).__init__(flowlines, mb_model=mb_model, y0=y0, inplace=inplace, - mb_elev_feedback='annual', **kwargs) - self.option_areaconstant = option_areaconstant - self.constantarea_years = constantarea_years - self.spinupyears = spinupyears - self.glac_idx_initial = [fl.thick.nonzero()[0] for fl in flowlines] - self.y0 = 0 - self.is_tidewater = is_tidewater - self.water_level = water_level + 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. -# widths_t0 = flowlines[0].widths_m -# area_v1 = widths_t0 * flowlines[0].dx_meter -# print('area v1:', area_v1.sum()) -# area_v2 = np.copy(area_v1) -# area_v2[flowlines[0].thick == 0] = 0 -# print('area v2:', area_v2.sum()) - - # HERE IS THE STUFF TO RECORD FOR EACH FLOWLINE! - if self.is_tidewater: - self.calving_k = cfg.PARAMS['calving_k'] - self.calving_m3_since_y0 = 0. # total calving since time y0 - - assert len(flowlines) == 1, 'MassRedistributionCurveModel is not set up for multiple flowlines' - - - def run_until(self, y1, run_single_year=False): - """Runs the model from the current year up to a given year date y1. + 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 - This function runs the model for the time difference y1-self.y0 - If self.y0 has not been specified at some point, it is 0 and y1 will - be the time span in years to run the model for. + 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 - Parameters - ---------- - y1 : float - Upper time span for how long the model should run - """ - - # We force timesteps to yearly timesteps - if run_single_year: - self.updategeometry(y1) - else: - years = np.arange(self.yr, y1) - for year in years: - self.updategeometry(year) - - # Check for domain bounds - if self.check_for_boundaries: - if self.fls[-1].thick[-1] > 10: - raise RuntimeError('Glacier exceeds domain boundaries, ' - 'at year: {}'.format(self.yr)) - # Check for NaNs - for fl in self.fls: - if np.any(~np.isfinite(fl.thick)): - raise FloatingPointError('NaN in numerical solution.') - - + 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 - def run_until_and_store(self, y1, run_path=None, diag_path=None, - store_monthly_step=None): - """Runs the model and returns intermediate steps in xarray datasets. + 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 - This function repeatedly calls FlowlineModel.run_until for either - monthly or yearly time steps up till the upper time boundary y1. + # --------------------------------------------------------------------------- + # 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 ---------- - y1 : int - Upper time span for how long the model should run (needs to be - a full year) - run_path : str - Path and filename where to store the model run dataset - diag_path : str - Path and filename where to store the model diagnostics dataset - store_monthly_step : Bool - If True (False) model diagnostics will be stored monthly (yearly). - If unspecified, we follow the update of the MB model, which - defaults to yearly (see __init__). + 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 ------- - run_ds : xarray.Dataset - stores the entire glacier geometry. It is useful to visualize the - glacier geometry or to restart a new run from a modelled geometry. - The glacier state is stored at the begining of each hydrological - year (not in between in order to spare disk space). - diag_ds : xarray.Dataset - stores a few diagnostic variables such as the volume, area, length - and ELA of the glacier. + 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 - if int(y1) != y1: - raise InvalidParamsError('run_until_and_store only accepts ' - 'integer year dates.') - - if not self.mb_model.hemisphere: - raise InvalidParamsError('run_until_and_store needs a ' - 'mass-balance model with an unambiguous ' - 'hemisphere.') - # time - yearly_time = np.arange(np.floor(self.yr), np.floor(y1)+1) - if store_monthly_step is None: - store_monthly_step = self.mb_step == 'monthly' + # --------------------------------------------------------------------------- + # calculate the yield ice thickness - if store_monthly_step: - monthly_time = utils.monthly_timeseries(self.yr, y1) - else: - monthly_time = np.arange(np.floor(self.yr), np.floor(y1)+1) + def balance_thickness(yield_strength, bed_elev): + """ + Ice thickness such that the stress matches the yield strength. - sm = cfg.PARAMS['hydro_month_' + self.mb_model.hemisphere] + 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 - yrs, months = utils.floatyear_to_date(monthly_time) - cyrs, cmonths = utils.hydrodate_to_calendardate(yrs, months, - start_month=sm) + 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)) ** 1)) + # TODO: Check on exponent on last term. In Ultee & Bassis 2016, this is squared, but in Ultee & Bassis 2020 supplement, it isn't. - # init output - if run_path is not None: - self.to_netcdf(run_path) - ny = len(yearly_time) - if ny == 1: - yrs = [yrs] - cyrs = [cyrs] - months = [months] - cmonths = [cmonths] - nm = len(monthly_time) - sects = [(np.zeros((ny, fl.nx)) * np.NaN) for fl in self.fls] - widths = [(np.zeros((ny, fl.nx)) * np.NaN) for fl in self.fls] - bucket = [(np.zeros(ny) * np.NaN) for _ in self.fls] - diag_ds = xr.Dataset() + # --------------------------------------------------------------------------- + # 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 - # Global attributes - diag_ds.attrs['description'] = 'OGGM model output' - diag_ds.attrs['oggm_version'] = __version__ - diag_ds.attrs['calendar'] = '365-day no leap' - diag_ds.attrs['creation_date'] = strftime("%Y-%m-%d %H:%M:%S", - gmtime()) - diag_ds.attrs['hemisphere'] = self.mb_model.hemisphere - diag_ds.attrs['water_level'] = self.water_level + # 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) - # Coordinates - diag_ds.coords['time'] = ('time', monthly_time) - diag_ds.coords['hydro_year'] = ('time', yrs) - diag_ds.coords['hydro_month'] = ('time', months) - diag_ds.coords['calendar_year'] = ('time', cyrs) - diag_ds.coords['calendar_month'] = ('time', cmonths) + 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 ") + - diag_ds['time'].attrs['description'] = 'Floating hydrological year' - diag_ds['hydro_year'].attrs['description'] = 'Hydrological year' - diag_ds['hydro_month'].attrs['description'] = 'Hydrological month' - diag_ds['calendar_year'].attrs['description'] = 'Calendar year' - diag_ds['calendar_month'].attrs['description'] = 'Calendar month' + 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 - # Variables and attributes - diag_ds['volume_m3'] = ('time', np.zeros(nm) * np.NaN) - diag_ds['volume_m3'].attrs['description'] = 'Total glacier volume' - diag_ds['volume_m3'].attrs['unit'] = 'm 3' - if self.is_tidewater: - diag_ds['volume_bsl_m3'] = ('time', np.zeros(nm) * np.NaN) - diag_ds['volume_bsl_m3'].attrs['description'] = ('Glacier volume ' - 'below ' - 'sea-level') - diag_ds['volume_bsl_m3'].attrs['unit'] = 'm 3' - diag_ds['volume_bwl_m3'] = ('time', np.zeros(nm) * np.NaN) - diag_ds['volume_bwl_m3'].attrs['description'] = ('Glacier volume ' - 'below ') - diag_ds['volume_bwl_m3'].attrs['unit'] = 'm 3' - diag_ds['area_m2'] = ('time', np.zeros(nm) * np.NaN) - diag_ds['area_m2'].attrs['description'] = 'Total glacier area' - diag_ds['area_m2'].attrs['unit'] = 'm 2' - diag_ds['length_m'] = ('time', np.zeros(nm) * np.NaN) - diag_ds['length_m'].attrs['description'] = 'Glacier length' - diag_ds['length_m'].attrs['unit'] = 'm 3' - diag_ds['ela_m'] = ('time', np.zeros(nm) * np.NaN) - diag_ds['ela_m'].attrs['description'] = ('Annual Equilibrium Line ' - 'Altitude (ELA)') - diag_ds['ela_m'].attrs['unit'] = 'm a.s.l' - if self.is_tidewater: - diag_ds['calving_m3'] = ('time', np.zeros(nm) * np.NaN) - diag_ds['calving_m3'].attrs['description'] = ('Total accumulated ' - 'calving flux') - diag_ds['calving_m3'].attrs['unit'] = 'm 3' - diag_ds['calving_rate_myr'] = ('time', np.zeros(nm) * np.NaN) - diag_ds['calving_rate_myr'].attrs['description'] = 'Calving rate' - diag_ds['calving_rate_myr'].attrs['unit'] = 'm yr-1' - # Run - j = 0 - for i, (yr, mo) in enumerate(zip(yearly_time[:-1], months[:-1])): - - # Record initial parameters - if i == 0: - diag_ds['volume_m3'].data[i] = self.volume_m3 - diag_ds['area_m2'].data[i] = self.area_m2 - diag_ds['length_m'].data[i] = self.length_m - - if self.is_tidewater: - diag_ds['volume_bsl_m3'].data[i] = self.volume_bsl_m3 - diag_ds['volume_bwl_m3'].data[i] = self.volume_bwl_m3 - - self.run_until(yr, run_single_year=True) - # Model run - if mo == 1: - for s, w, b, fl in zip(sects, widths, bucket, self.fls): - s[j, :] = fl.section - w[j, :] = fl.widths_m - if self.is_tidewater: - try: - b[j] = fl.calving_bucket_m3 - except AttributeError: - pass - j += 1 - # Diagnostics - diag_ds['volume_m3'].data[i+1] = self.volume_m3 - diag_ds['area_m2'].data[i+1] = self.area_m2 - diag_ds['length_m'].data[i+1] = self.length_m - if self.is_tidewater: - diag_ds['calving_m3'].data[i+1] = self.calving_m3_since_y0 - diag_ds['calving_rate_myr'].data[i+1] = self.calving_rate_myr - diag_ds['volume_bsl_m3'].data[i+1] = self.volume_bsl_m3 - diag_ds['volume_bwl_m3'].data[i+1] = self.volume_bwl_m3 - # to datasets - run_ds = [] - for (s, w, b) in zip(sects, widths, bucket): - ds = xr.Dataset() - ds.attrs['description'] = 'OGGM model output' - ds.attrs['oggm_version'] = __version__ - ds.attrs['calendar'] = '365-day no leap' - ds.attrs['creation_date'] = strftime("%Y-%m-%d %H:%M:%S", - gmtime()) - ds.coords['time'] = yearly_time - ds['time'].attrs['description'] = 'Floating hydrological year' - varcoords = OrderedDict(time=('time', yearly_time), - year=('time', yearly_time)) - ds['ts_section'] = xr.DataArray(s, dims=('time', 'x'), - coords=varcoords) - ds['ts_width_m'] = xr.DataArray(w, dims=('time', 'x'), - coords=varcoords) - if self.is_tidewater: - ds['ts_calving_bucket_m3'] = xr.DataArray(b, dims=('time', ), - coords=varcoords) - run_ds.append(ds) - # write output? - if run_path is not None: - encode = {'ts_section': {'zlib': True, 'complevel': 5}, - 'ts_width_m': {'zlib': True, 'complevel': 5}, - } - for i, ds in enumerate(run_ds): - ds.to_netcdf(run_path, 'a', group='fl_{}'.format(i), - encoding=encode) - # Add other diagnostics - diag_ds.to_netcdf(run_path, 'a') - if diag_path is not None: - diag_ds.to_netcdf(diag_path) - return run_ds, diag_ds - - - def updategeometry(self, year, debug=False): - """Update geometry for a given year""" +#%% +class MassRedistributionCurveModel(FlowlineModel): + """Glacier geometry updated using mass redistribution curves; also known as the "delta-h method" + + This uses mass redistribution curves from Huss et al. (2010) to update the glacier geometry + """ + + def __init__(self, flowlines, mb_model=None, y0=0., + glen_a=None, fs=0., is_tidewater=False, water_level=0, +# calving_k=0, + inplace=False, + debug=True, + option_areaconstant=False, spinupyears=0, + constantarea_years=0, + **kwargs): + """ Instanciate the model. - if debug: - print('year:', year) - - # Loop over flowlines - for fl_id, fl in enumerate(self.fls): + Parameters + ---------- + flowlines : list + the glacier flowlines + mb_model : MassBalanceModel + the mass-balance model + y0 : int + initial year of the simulation + inplace : bool + whether or not to make a copy of the flowline objects for the run + setting to True implies that your objects will be modified at run + time by the model (can help to spare memory) + is_tidewater: bool, default: False + use the very basic parameterization for tidewater glaciers + mb_elev_feedback : str, default: 'annual' + 'never', 'always', 'annual', or 'monthly': how often the + mass-balance should be recomputed from the mass balance model. + 'Never' is equivalent to 'annual' but without elevation feedback + at all (the heights are taken from the first call). + check_for_boundaries: bool, default: True + raise an error when the glacier grows bigger than the domain + boundaries + """ + super(MassRedistributionCurveModel, self).__init__(flowlines, mb_model=mb_model, y0=y0, inplace=inplace, + mb_elev_feedback='annual', **kwargs) + self.option_areaconstant = option_areaconstant + self.constantarea_years = constantarea_years + self.spinupyears = spinupyears + self.glac_idx_initial = [fl.thick.nonzero()[0] for fl in flowlines] + self.y0 = 0 + self.is_tidewater = is_tidewater + self.water_level = water_level - # Flowline state - heights = self.fls[fl_id].surface_h.copy() - section_t0 = self.fls[fl_id].section.copy() - thick_t0 = self.fls[fl_id].thick.copy() - width_t0 = self.fls[fl_id].widths_m.copy() - - # CONSTANT AREAS - # Mass redistribution ignored for calibration and spinup years (glacier properties constant) - if (self.option_areaconstant) or (year < self.spinupyears) or (year < self.constantarea_years): - # run mass balance - glac_bin_massbalclim_annual = self.mb_model.get_annual_mb(heights, fls=self.fls, fl_id=fl_id, - year=year, debug=False) - # MASS REDISTRIBUTION - else: - # FRONTAL ABLATION - if self.is_tidewater: - # Frontal ablation (m3 ice) - fa_m3 = self._get_annual_frontalablation(heights, fls=self.fls, fl_id=fl_id, - year=year, debug=False) - if debug: - print('fa_m3_init:', fa_m3) - vol_init = (self.fls[fl_id].section * fl.dx_meter).sum() - print(' volume init:', np.round(vol_init)) - print(' volume final:', np.round(vol_init-fa_m3)) - # First, remove volume lost to frontal ablation - # changes to _t0 not _t1, since t1 will be done in the mass redistribution - - glac_idx_bsl = np.where((thick_t0 > 0) & (fl.bed_h < self.water_level))[0] - while fa_m3 > 0 and len(glac_idx_bsl) > 0: - if debug: - print('fa_m3_remaining:', fa_m3) +# widths_t0 = flowlines[0].widths_m +# area_v1 = widths_t0 * flowlines[0].dx_meter +# print('area v1:', area_v1.sum()) +# area_v2 = np.copy(area_v1) +# area_v2[flowlines[0].thick == 0] = 0 +# print('area v2:', area_v2.sum()) + + # HERE IS THE STUFF TO RECORD FOR EACH FLOWLINE! + if self.is_tidewater: + self.calving_k = cfg.PARAMS['calving_k'] + self.calving_m3_since_y0 = 0. # total calving since time y0 + + assert len(flowlines) == 1, 'MassRedistributionCurveModel is not set up for multiple flowlines' + + + def run_until(self, y1, run_single_year=False): + """Runs the model from the current year up to a given year date y1. - # OGGM code -# glac_idx_bsl = np.where((thick_t0 > 0) & (fl.bed_h < self.water_level))[0] - last_idx = glac_idx_bsl[-1] - - if debug: - print('before:', np.round(self.fls[fl_id].section[last_idx],0), - np.round(self.fls[fl_id].thick[last_idx],0), - np.round(heights[last_idx],0)) - - vol_last = section_t0[last_idx] * fl.dx_meter - - # If frontal ablation more than bin volume, remove entire bin - if fa_m3 > vol_last: - # Record frontal ablation (m3 w.e.) in mass balance model for output - self.mb_model.glac_bin_frontalablation[last_idx,int(12*(year+1)-1)] = ( - vol_last * pygem_prms.density_ice / pygem_prms.density_water) - # Update ice thickness and section area - section_t0[last_idx] = 0 - self.fls[fl_id].section = section_t0 - # Remove volume from frontal ablation "bucket" - fa_m3 -= vol_last - - # Otherwise, remove ice from the section - else: - # Update section to remove frontal ablation - section_t0[last_idx] = section_t0[last_idx] - fa_m3 / fl.dx_meter - self.fls[fl_id].section = section_t0 - # Record frontal ablation(m3 w.e.) - self.mb_model.glac_bin_frontalablation[last_idx,int(12*(year+1)-1)] = ( - fa_m3 * pygem_prms.density_ice / pygem_prms.density_water) - # Frontal ablation bucket now empty - fa_m3 = 0 - - # Update flowline - heights = self.fls[fl_id].surface_h.copy() - section_t0 = self.fls[fl_id].section.copy() - thick_t0 = self.fls[fl_id].thick.copy() - width_t0 = self.fls[fl_id].widths_m.copy() - - if debug: - print('after:', np.round(self.fls[fl_id].section[last_idx],0), - np.round(self.fls[fl_id].thick[last_idx],0), - np.round(heights[last_idx],0)) - print(' vol final:', (self.fls[fl_id].section * fl.dx_meter).sum()) - - glac_idx_bsl = np.where((thick_t0 > 0) & (fl.bed_h < self.water_level))[0] - - - # Redistribute mass if glacier was not fully removed by frontal ablation - if len(section_t0.nonzero()[0]) > 0: - # Mass redistribution according to Huss empirical curves - # Annual glacier mass balance [m ice s-1] - glac_bin_massbalclim_annual = self.mb_model.get_annual_mb(heights, fls=self.fls, fl_id=fl_id, - year=year, debug=False) - sec_in_year = (self.mb_model.dates_table.loc[12*year:12*(year+1)-1,'daysinmonth'].values.sum() - * 24 * 3600) - -# print(' volume change [m3]:', (glac_bin_massbalclim_annual * sec_in_year * -# (width_t0 * fl.dx_meter)).sum()) -# print(glac_bin_masssbalclim_annual) -# print(sec_in_year) -# print(width_t0.sum()) -# print(fl.dx_meter) -# print(width_t0 * fl.dx_meter) - -# # Debugging block -# debug_years = [71] -# if year in debug_years: -# print(year, glac_bin_massbalclim_annual) -# print('section t0:', section_t0) -# print('thick_t0:', thick_t0) -# print('width_t0:', width_t0) -# print(self.glac_idx_initial[fl_id]) -# print('heights:', heights) + This function runs the model for the time difference y1-self.y0 + If self.y0 has not been specified at some point, it is 0 and y1 will + be the time span in years to run the model for. + + Parameters + ---------- + y1 : float + Upper time span for how long the model should run + """ - self._massredistributionHuss(section_t0, thick_t0, width_t0, glac_bin_massbalclim_annual, - self.glac_idx_initial[fl_id], heights, sec_in_year=sec_in_year) + # We force timesteps to yearly timesteps + if run_single_year: + self.updategeometry(y1) + else: + years = np.arange(self.yr, y1) + for year in years: + self.updategeometry(year) + + # Check for domain bounds + if self.check_for_boundaries: + if self.fls[-1].thick[-1] > 10: + raise RuntimeError('Glacier exceeds domain boundaries, ' + 'at year: {}'.format(self.yr)) + # Check for NaNs + for fl in self.fls: + if np.any(~np.isfinite(fl.thick)): + raise FloatingPointError('NaN in numerical solution.') + - # Record glacier properties (volume [m3], area [m2], thickness [m], width [km]) - # record the next year's properties as well - # 'year + 1' used so the glacier properties are consistent with mass balance computations - year = int(year) # required to ensure proper indexing with run_until_and_store (10/21/2020) - glacier_area = fl.widths_m * fl.dx_meter - glacier_area[fl.thick == 0] = 0 - self.mb_model.glac_bin_area_annual[:,year+1] = glacier_area - self.mb_model.glac_bin_icethickness_annual[:,year+1] = fl.thick - self.mb_model.glac_bin_width_annual[:,year+1] = fl.widths_m - self.mb_model.glac_wide_area_annual[year+1] = glacier_area.sum() - self.mb_model.glac_wide_volume_annual[year+1] = (fl.section * fl.dx_meter).sum() - - - #%% ----- Calving Law ---------- -def fa_sermeq_speed_law(model,last_above_wl, v_scaling=1, verbose=False, - tau0=150e3, 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 + def run_until_and_store(self, y1, run_path=None, diag_path=None, + store_monthly_step=None): + """Runs the model and returns intermediate steps in xarray datasets. - 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=150e3, 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. + This function repeatedly calls FlowlineModel.run_until for either + monthly or yearly time steps up till the upper time boundary y1. 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. + y1 : int + Upper time span for how long the model should run (needs to be + a full year) + run_path : str + Path and filename where to store the model run dataset + diag_path : str + Path and filename where to store the model diagnostics dataset + store_monthly_step : Bool + If True (False) model diagnostics will be stored monthly (yearly). + If unspecified, we follow the update of the MB model, which + defaults to yearly (see __init__). Returns ------- - tau_y: float - The yield strength for these conditions. + run_ds : xarray.Dataset + stores the entire glacier geometry. It is useful to visualize the + glacier geometry or to restart a new run from a modelled geometry. + The glacier state is stored at the begining of each hydrological + year (not in between in order to spare disk space). + diag_ds : xarray.Dataset + stores a few diagnostic variables such as the volume, area, length + and ELA of the glacier. """ - 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 = tau0 + mu * N - else: # assume constant if not set - ty = tau0 - return ty - - # --------------------------------------------------------------------------- - # calculate the yield ice thickness + if int(y1) != y1: + raise InvalidParamsError('run_until_and_store only accepts ' + 'integer year dates.') - def balance_thickness(yield_strength, bed_elev): - """ - Ice thickness such that the stress matches the yield strength. + if not self.mb_model.hemisphere: + raise InvalidParamsError('run_until_and_store needs a ' + 'mass-balance model with an unambiguous ' + 'hemisphere.') + # time + yearly_time = np.arange(np.floor(self.yr), np.floor(y1)+1) - 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 + if store_monthly_step is None: + store_monthly_step = self.mb_step == 'monthly' - Returns - ------- - Hy: float - The ice thickness for stress balance at the terminus. - """ - if bed_elev < 0: - D = -1 * bed_elev + if store_monthly_step: + monthly_time = utils.monthly_timeseries(self.yr, y1) else: - D = 0 - return (2 * yield_strength / (RHO_ICE * G)) + np.sqrt( - (RHO_SEA * (D ** 2) / RHO_ICE) + ((2 * yield_strength / (RHO_ICE * G)) ** 1)) - # 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 + monthly_time = np.arange(np.floor(self.yr), np.floor(y1)+1) + + sm = cfg.PARAMS['hydro_month_' + self.mb_model.hemisphere] + + yrs, months = utils.floatyear_to_date(monthly_time) + cyrs, cmonths = utils.hydrodate_to_calendardate(yrs, months, + start_month=sm) + + # init output + if run_path is not None: + self.to_netcdf(run_path) + ny = len(yearly_time) + if ny == 1: + yrs = [yrs] + cyrs = [cyrs] + months = [months] + cmonths = [cmonths] + nm = len(monthly_time) + sects = [(np.zeros((ny, fl.nx)) * np.NaN) for fl in self.fls] + widths = [(np.zeros((ny, fl.nx)) * np.NaN) for fl in self.fls] + bucket = [(np.zeros(ny) * np.NaN) for _ in self.fls] + diag_ds = xr.Dataset() + + # Global attributes + diag_ds.attrs['description'] = 'OGGM model output' + diag_ds.attrs['oggm_version'] = __version__ + diag_ds.attrs['calendar'] = '365-day no leap' + diag_ds.attrs['creation_date'] = strftime("%Y-%m-%d %H:%M:%S", + gmtime()) + diag_ds.attrs['hemisphere'] = self.mb_model.hemisphere + diag_ds.attrs['water_level'] = self.water_level + + # Coordinates + diag_ds.coords['time'] = ('time', monthly_time) + diag_ds.coords['hydro_year'] = ('time', yrs) + diag_ds.coords['hydro_month'] = ('time', months) + diag_ds.coords['calendar_year'] = ('time', cyrs) + diag_ds.coords['calendar_month'] = ('time', cmonths) + + diag_ds['time'].attrs['description'] = 'Floating hydrological year' + diag_ds['hydro_year'].attrs['description'] = 'Hydrological year' + diag_ds['hydro_month'].attrs['description'] = 'Hydrological month' + diag_ds['calendar_year'].attrs['description'] = 'Calendar year' + diag_ds['calendar_month'].attrs['description'] = 'Calendar month' + + # Variables and attributes + diag_ds['volume_m3'] = ('time', np.zeros(nm) * np.NaN) + diag_ds['volume_m3'].attrs['description'] = 'Total glacier volume' + diag_ds['volume_m3'].attrs['unit'] = 'm 3' + if self.is_tidewater: + diag_ds['volume_bsl_m3'] = ('time', np.zeros(nm) * np.NaN) + diag_ds['volume_bsl_m3'].attrs['description'] = ('Glacier volume ' + 'below ' + 'sea-level') + diag_ds['volume_bsl_m3'].attrs['unit'] = 'm 3' + diag_ds['volume_bwl_m3'] = ('time', np.zeros(nm) * np.NaN) + diag_ds['volume_bwl_m3'].attrs['description'] = ('Glacier volume ' + 'below ') + diag_ds['volume_bwl_m3'].attrs['unit'] = 'm 3' + + diag_ds['area_m2'] = ('time', np.zeros(nm) * np.NaN) + diag_ds['area_m2'].attrs['description'] = 'Total glacier area' + diag_ds['area_m2'].attrs['unit'] = 'm 2' + diag_ds['length_m'] = ('time', np.zeros(nm) * np.NaN) + diag_ds['length_m'].attrs['description'] = 'Glacier length' + diag_ds['length_m'].attrs['unit'] = 'm 3' + diag_ds['ela_m'] = ('time', np.zeros(nm) * np.NaN) + diag_ds['ela_m'].attrs['description'] = ('Annual Equilibrium Line ' + 'Altitude (ELA)') + diag_ds['ela_m'].attrs['unit'] = 'm a.s.l' + if self.is_tidewater: + diag_ds['calving_m3'] = ('time', np.zeros(nm) * np.NaN) + diag_ds['calving_m3'].attrs['description'] = ('Total accumulated ' + 'calving flux') + diag_ds['calving_m3'].attrs['unit'] = 'm 3' + diag_ds['calving_rate_myr'] = ('time', np.zeros(nm) * np.NaN) + diag_ds['calving_rate_myr'].attrs['description'] = 'Calving rate' + diag_ds['calving_rate_myr'].attrs['unit'] = 'm yr-1' + + # Run + j = 0 + for i, (yr, mo) in enumerate(zip(yearly_time[:-1], months[:-1])): + + # Record initial parameters + if i == 0: + diag_ds['volume_m3'].data[i] = self.volume_m3 + diag_ds['area_m2'].data[i] = self.area_m2 + diag_ds['length_m'].data[i] = self.length_m + + if self.is_tidewater: + diag_ds['volume_bsl_m3'].data[i] = self.volume_bsl_m3 + diag_ds['volume_bwl_m3'].data[i] = self.volume_bwl_m3 + + self.run_until(yr, run_single_year=True) + # Model run + if mo == 1: + for s, w, b, fl in zip(sects, widths, bucket, self.fls): + s[j, :] = fl.section + w[j, :] = fl.widths_m + if self.is_tidewater: + try: + b[j] = fl.calving_bucket_m3 + except AttributeError: + pass + j += 1 + # Diagnostics + diag_ds['volume_m3'].data[i+1] = self.volume_m3 + diag_ds['area_m2'].data[i+1] = self.area_m2 + diag_ds['length_m'].data[i+1] = self.length_m + + if self.is_tidewater: + diag_ds['calving_m3'].data[i+1] = self.calving_m3_since_y0 + diag_ds['calving_rate_myr'].data[i+1] = self.calving_rate_myr + diag_ds['volume_bsl_m3'].data[i+1] = self.volume_bsl_m3 + diag_ds['volume_bwl_m3'].data[i+1] = self.volume_bwl_m3 + + # to datasets + run_ds = [] + for (s, w, b) in zip(sects, widths, bucket): + ds = xr.Dataset() + ds.attrs['description'] = 'OGGM model output' + ds.attrs['oggm_version'] = __version__ + ds.attrs['calendar'] = '365-day no leap' + ds.attrs['creation_date'] = strftime("%Y-%m-%d %H:%M:%S", + gmtime()) + ds.coords['time'] = yearly_time + ds['time'].attrs['description'] = 'Floating hydrological year' + varcoords = OrderedDict(time=('time', yearly_time), + year=('time', yearly_time)) + ds['ts_section'] = xr.DataArray(s, dims=('time', 'x'), + coords=varcoords) + ds['ts_width_m'] = xr.DataArray(w, dims=('time', 'x'), + coords=varcoords) + if self.is_tidewater: + ds['ts_calving_bucket_m3'] = xr.DataArray(b, dims=('time', ), + coords=varcoords) + run_ds.append(ds) + + # write output? + if run_path is not None: + encode = {'ts_section': {'zlib': True, 'complevel': 5}, + 'ts_width_m': {'zlib': True, 'complevel': 5}, + } + for i, ds in enumerate(run_ds): + ds.to_netcdf(run_path, 'a', group='fl_{}'.format(i), + encoding=encode) + # Add other diagnostics + diag_ds.to_netcdf(run_path, 'a') + + if diag_path is not None: + diag_ds.to_netcdf(diag_path) + + return run_ds, diag_ds + + + def updategeometry(self, year, debug=False): + """Update geometry for a given year""" + + if debug: + print('year:', year) + + # Loop over flowlines + for fl_id, fl in enumerate(self.fls): + + # Flowline state + heights = self.fls[fl_id].surface_h.copy() + section_t0 = self.fls[fl_id].section.copy() + thick_t0 = self.fls[fl_id].thick.copy() + width_t0 = self.fls[fl_id].widths_m.copy() + + # CONSTANT AREAS + # Mass redistribution ignored for calibration and spinup years (glacier properties constant) + if (self.option_areaconstant) or (year < self.spinupyears) or (year < self.constantarea_years): + # run mass balance + glac_bin_massbalclim_annual = self.mb_model.get_annual_mb(heights, fls=self.fls, fl_id=fl_id, + year=year, debug=False) + # MASS REDISTRIBUTION + else: + # FRONTAL ABLATION + if self.is_tidewater: + # Frontal ablation (m3 ice) + fa_m3 = self._get_annual_frontalablation(heights, fls=self.fls, fl_id=fl_id, + year=year, debug=False) + if debug: + print('fa_m3_init:', fa_m3) + vol_init = (self.fls[fl_id].section * fl.dx_meter).sum() + print(' volume init:', np.round(vol_init)) + print(' volume final:', np.round(vol_init-fa_m3)) + # First, remove volume lost to frontal ablation + # changes to _t0 not _t1, since t1 will be done in the mass redistribution + + glac_idx_bsl = np.where((thick_t0 > 0) & (fl.bed_h < self.water_level))[0] + while fa_m3 > 0 and len(glac_idx_bsl) > 0: + if debug: + print('fa_m3_remaining:', fa_m3) + + # OGGM code +# glac_idx_bsl = np.where((thick_t0 > 0) & (fl.bed_h < self.water_level))[0] + last_idx = glac_idx_bsl[-1] + + if debug: + print('before:', np.round(self.fls[fl_id].section[last_idx],0), + np.round(self.fls[fl_id].thick[last_idx],0), + np.round(heights[last_idx],0)) + + vol_last = section_t0[last_idx] * fl.dx_meter + + # If frontal ablation more than bin volume, remove entire bin + if fa_m3 > vol_last: + # Record frontal ablation (m3 w.e.) in mass balance model for output + self.mb_model.glac_bin_frontalablation[last_idx,int(12*(year+1)-1)] = ( + vol_last * pygem_prms.density_ice / pygem_prms.density_water) + # Update ice thickness and section area + section_t0[last_idx] = 0 + self.fls[fl_id].section = section_t0 + # Remove volume from frontal ablation "bucket" + fa_m3 -= vol_last + + # Otherwise, remove ice from the section + else: + # Update section to remove frontal ablation + section_t0[last_idx] = section_t0[last_idx] - fa_m3 / fl.dx_meter + self.fls[fl_id].section = section_t0 + # Record frontal ablation(m3 w.e.) + self.mb_model.glac_bin_frontalablation[last_idx,int(12*(year+1)-1)] = ( + fa_m3 * pygem_prms.density_ice / pygem_prms.density_water) + # Frontal ablation bucket now empty + fa_m3 = 0 + + # Update flowline + heights = self.fls[fl_id].surface_h.copy() + section_t0 = self.fls[fl_id].section.copy() + thick_t0 = self.fls[fl_id].thick.copy() + width_t0 = self.fls[fl_id].widths_m.copy() + + if debug: + print('after:', np.round(self.fls[fl_id].section[last_idx],0), + np.round(self.fls[fl_id].thick[last_idx],0), + np.round(heights[last_idx],0)) + print(' vol final:', (self.fls[fl_id].section * fl.dx_meter).sum()) + + glac_idx_bsl = np.where((thick_t0 > 0) & (fl.bed_h < self.water_level))[0] + + + # Redistribute mass if glacier was not fully removed by frontal ablation + if len(section_t0.nonzero()[0]) > 0: + # Mass redistribution according to Huss empirical curves + # Annual glacier mass balance [m ice s-1] + glac_bin_massbalclim_annual = self.mb_model.get_annual_mb(heights, fls=self.fls, fl_id=fl_id, + year=year, debug=False) + sec_in_year = (self.mb_model.dates_table.loc[12*year:12*(year+1)-1,'daysinmonth'].values.sum() + * 24 * 3600) + +# print(' volume change [m3]:', (glac_bin_massbalclim_annual * sec_in_year * +# (width_t0 * fl.dx_meter)).sum()) +# print(glac_bin_masssbalclim_annual) +# print(sec_in_year) +# print(width_t0.sum()) +# print(fl.dx_meter) +# print(width_t0 * fl.dx_meter) + +# # Debugging block +# debug_years = [71] +# if year in debug_years: +# print(year, glac_bin_massbalclim_annual) +# print('section t0:', section_t0) +# print('thick_t0:', thick_t0) +# print('width_t0:', width_t0) +# print(self.glac_idx_initial[fl_id]) +# print('heights:', heights) + + self._massredistributionHuss(section_t0, thick_t0, width_t0, glac_bin_massbalclim_annual, + self.glac_idx_initial[fl_id], heights, sec_in_year=sec_in_year) + + # Record glacier properties (volume [m3], area [m2], thickness [m], width [km]) + # record the next year's properties as well + # 'year + 1' used so the glacier properties are consistent with mass balance computations + year = int(year) # required to ensure proper indexing with run_until_and_store (10/21/2020) + glacier_area = fl.widths_m * fl.dx_meter + glacier_area[fl.thick == 0] = 0 + self.mb_model.glac_bin_area_annual[:,year+1] = glacier_area + self.mb_model.glac_bin_icethickness_annual[:,year+1] = fl.thick + self.mb_model.glac_bin_width_annual[:,year+1] = fl.widths_m + self.mb_model.glac_wide_area_annual[year+1] = glacier_area.sum() + self.mb_model.glac_wide_volume_annual[year+1] = (fl.section * fl.dx_meter).sum() - # 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 #%% ----- FRONTAL ABLATION ----- @@ -812,7 +823,7 @@ def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, #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,v_scaling = 1, verbose = False,tau0 = self.calving_k, + 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********************') diff --git a/pygem/pygem_modelsetup.py b/pygem/pygem_modelsetup.py index c09339b9..50ecff90 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -297,9 +297,12 @@ 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] @@ -308,6 +311,7 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' 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: @@ -318,8 +322,10 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' 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)) @@ -347,6 +353,7 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' axis=0)) 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 +400,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 From ec352b7785730b4795ef049cf5f99510d8788365 Mon Sep 17 00:00:00 2001 From: ruitang yang Date: Wed, 13 Dec 2023 15:48:27 +0100 Subject: [PATCH 06/45] modify the function fa_sermeq_speed_law --- pygem/glacierdynamics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index 6d583743..c79132e1 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -167,7 +167,7 @@ def balance_thickness(yield_strength, 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)) ** 1)) + (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. # --------------------------------------------------------------------------- From 1f544d7cd527365af12aa47389103f9d3a5a7f2a Mon Sep 17 00:00:00 2001 From: ruitang yang Date: Fri, 22 Dec 2023 14:42:32 +0100 Subject: [PATCH 07/45] print the filepath of mb_obs --- pygem/shop/mbdata.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pygem/shop/mbdata.py b/pygem/shop/mbdata.py index 39e3093d..84b51a10 100755 --- a/pygem/shop/mbdata.py +++ b/pygem/shop/mbdata.py @@ -111,6 +111,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) From d8ff180b79832fb103fd327a40fa961e2cc88c9c Mon Sep 17 00:00:00 2001 From: ruitang yang Date: Fri, 29 Dec 2023 16:44:46 +0100 Subject: [PATCH 08/45] add debug print in the function get_annual_mb --- pygem/massbalance.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 3cb9347a..5e9c3752 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -215,7 +215,12 @@ def get_annual_mb(self, heights, year=None, fls=None, fl_id=None, 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 @@ -225,10 +230,11 @@ def get_annual_mb(self, heights, year=None, fls=None, fl_id=None, # 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)) @@ -251,6 +257,7 @@ def get_annual_mb(self, heights, year=None, fls=None, fl_id=None, 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] + print('****************** off glacier idx is*******************',offglac_idx) # Functions currently set up for monthly timestep # only compute mass balance while glacier exists From fa77525157204593451ea415ee412a9f1827e6d9 Mon Sep 17 00:00:00 2001 From: ruitang yang Date: Tue, 16 Jan 2024 19:47:48 +0100 Subject: [PATCH 09/45] add the get_monthly_mb function; introduce new parameters in year_month in get_annual_mb; add the volume_m3_annual/month_ice in the function ensure_mass_conservation --- pygem/massbalance.py | 42 +++++++++++++++++++++++++++++++++++------- 1 file changed, 35 insertions(+), 7 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 5e9c3752..b0784134 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -180,9 +180,18 @@ 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_monthly_mb(self, heights, year=None, fls=None, fl_id=None, + debug=False, option_areaconstant=False): + year_floor=np.floor(year) + month=year-year_floor + + 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=None, - debug=False, option_areaconstant=False): + debug=False, option_areaconstant=False, year_month=None): """FIXED FORMAT FOR THE FLOWLINE MODEL Returns annual climatic mass balance [m ice per second] @@ -198,8 +207,9 @@ 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) @@ -631,10 +641,20 @@ def get_annual_mb(self, heights, year=None, fls=None, fl_id=None, ## 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 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: + seconds_in_month = self.dayspermonth[12*year+int(year_month*12)]* 24 * 3600 + 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(12*year+int(year_month*12)) + print(year_month) + print(year) + if self.inversion_filter: mb = np.minimum.accumulate(mb) @@ -837,7 +857,15 @@ def ensure_mass_conservation(self, diag): 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] + # if the dynamic step is monthly, the volume is monthly, need to be calculated as annual + Dynamic_step_Monthly =True + if Dynamic_step_Monthly : + calving_m3_month_ice = diag.calving_m3.values[1:] - diag.calving_m3.values[0:-1] + calving_m3_annual_ice = calving_m3_month_ice.reshape(-1,12).sum(1) + vol_change_annual_diag = calving_m3_annual_ice + 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 From 2110d67caaf62236d7b6aca50818615c59c89bb3 Mon Sep 17 00:00:00 2001 From: ruitang yang Date: Tue, 16 Jan 2024 19:50:19 +0100 Subject: [PATCH 10/45] modify the index for mb and volume_m3_month/annual_ice --- pygem/massbalance.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index b0784134..f135e0ac 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -182,8 +182,9 @@ def __init__(self, gdir, modelprms, glacier_rgi_table, 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-year_floor + month=year-int(year_floor) mb=self.get_annual_mb(heights=heights, year=year_floor, fls=fls, fl_id=fl_id, debug=debug, option_areaconstant=False, year_month=month) @@ -646,8 +647,8 @@ def get_annual_mb(self, heights, year=None, fls=None, fl_id=None, mb = (self.glac_bin_massbalclim[:,12*year:12*(year+1)].sum(1) * pygem_prms.density_water / pygem_prms.density_ice / seconds_in_year) else: - seconds_in_month = self.dayspermonth[12*year+int(year_month*12)]* 24 * 3600 - mb = (self.glac_bin_massbalclim[:,12*year+int(year_month*12)] + seconds_in_month = self.dayspermonth[12*year+round(year_month*12)]* 24 * 3600 + mb = (self.glac_bin_massbalclim[:,12*year+round(year_month*12)] * pygem_prms.density_water / pygem_prms.density_ice /seconds_in_month) print("index for mb") @@ -857,15 +858,18 @@ def ensure_mass_conservation(self, diag): 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) + #%% Dynamic running step # if the dynamic step is monthly, the volume is monthly, need to be calculated as annual Dynamic_step_Monthly =True if Dynamic_step_Monthly : - calving_m3_month_ice = diag.calving_m3.values[1:] - diag.calving_m3.values[0:-1] - calving_m3_annual_ice = calving_m3_month_ice.reshape(-1,12).sum(1) - vol_change_annual_diag = calving_m3_annual_ice + 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 From 1e954dd28dc638cbb1f2b67c96b9c07bbde2a9a6 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Sat, 20 Jan 2024 00:00:07 +0100 Subject: [PATCH 11/45] add variables length change and length annual in PyGEMMassBalance, but actually it's the wrong place, should be in OGGM --- pygem/massbalance.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index f135e0ac..bebbdfd2 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -120,6 +120,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,6 +139,7 @@ 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_wide_length_annual = np.zeros(self.nyears+1) self.glac_wide_massbaltotal = np.zeros(self.nmonths) self.glac_wide_runoff = np.zeros(self.nmonths) self.glac_wide_snowline = np.zeros(self.nmonths) From 93a11411c9d039fbf94e993c5a6a2819a038111c Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Sat, 20 Jan 2024 00:01:00 +0100 Subject: [PATCH 12/45] comment the length change and length annual variables --- pygem/massbalance.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index bebbdfd2..fb6e98f4 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -120,7 +120,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_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)) @@ -139,7 +139,7 @@ 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_wide_length_annual = np.zeros(self.nyears+1) + #self.glac_wide_length_annual = np.zeros(self.nyears+1) self.glac_wide_massbaltotal = np.zeros(self.nmonths) self.glac_wide_runoff = np.zeros(self.nmonths) self.glac_wide_snowline = np.zeros(self.nmonths) From 2229aea648970c1ac0c4cb89855460ae5c603be9 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Tue, 30 Jan 2024 21:58:50 +0100 Subject: [PATCH 13/45] add the debug input in PyGEMMassBalance --- pygem/massbalance.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index fb6e98f4..190c7055 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -48,6 +48,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("PYGEM: flowlines") + print(type(fls).__name__) + if debug: print('\n\nDEBUGGING MASS BALANCE FUNCTION\n\n') self.debug_refreeze = debug_refreeze From a7996cef38f66ff073a144b2a7e26de470c5908f Mon Sep 17 00:00:00 2001 From: "Brandon S. Tober" Date: Wed, 3 Apr 2024 10:33:37 -0400 Subject: [PATCH 14/45] Output class objects (#39) * bug fix in pygem_modelsetup.datesmodelrun(). pandas.date_range() has a timestamp range limitation based on the chosen resolution. The deafult is nanoseconds, which was reduced to seconds: https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#timeseries-timestamp-limits * check for differing climate data file paths or names based * framework for PyGEM output classes * added conda env file * output dataclasses * output cleanup * output binned stats * model setup group by thousands * output classes progress * Delete environment.yml --------- Co-authored-by: David Rounce --- pygem/output.py | 651 ++++++++++++++++++++++++++++++++++++++ pygem/pygem_modelsetup.py | 13 +- 2 files changed, 659 insertions(+), 5 deletions(-) create mode 100644 pygem/output.py diff --git a/pygem/output.py b/pygem/output.py new file mode 100644 index 00000000..1decf628 --- /dev/null +++ b/pygem/output.py @@ -0,0 +1,651 @@ +#!/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 cftime +import os +import 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 + wateryear : bool + pygem_version : float + user_info : dict + outdir : str + gcm_name : str + scenario : str + realization : str + calib_opt : int + sim_iters : int + modelprms : dict + ba_opt : int + gcm_bc_startyr : int + gcm_startyr : int + gcm_endyr: 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.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 self.calib_opt: + self.outfn += f'{self.calib_opt}_' + 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{self.ba_opt}_' + yr0 = self.gcm_bc_startyr + else: + self.outfn += 'ba0_' + yr0 = self.gcm_startyr + if self.calib_opt: + self.outfn += 'SETS_' + self.outfn += f'{yr0}_' + self.outfn += f'{self.gcm_endyr}_' + + def get_fn(self): + return self.outfn + + def set_time_vals(self): + if self.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 self.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 self.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): + pass + + 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': self.user_info['institution'], + 'history': f'Created by {self.user_info["name"]} ({self.user_info["email"]}) on ' + datetime.today().strftime('%Y-%m-%d'), + 'references': 'doi:10.3389/feart.2019.00331 and doi:10.1017/jog.2019.91'} + + 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 + """ + extra_vars : bool + + 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 self.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..7c7a9253 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -264,7 +264,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. @@ -315,14 +316,16 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' # 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 else: glacier_table = pd.concat([glacier_table, csv_regionO1], axis=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] @@ -494,7 +497,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() From 79fa1f8bb32c40749b955a1ce5cf245786c14139 Mon Sep 17 00:00:00 2001 From: "Brandon S. Tober" Date: Wed, 3 Apr 2024 20:40:54 -0400 Subject: [PATCH 15/45] output cleanup, minor bug fixes (#40) --- pygem/output.py | 65 ++++++++++++++++++++++++++++--------------------- 1 file changed, 37 insertions(+), 28 deletions(-) diff --git a/pygem/output.py b/pygem/output.py index 1decf628..c499957f 100644 --- a/pygem/output.py +++ b/pygem/output.py @@ -19,9 +19,7 @@ import numpy as np import pandas as pd import xarray as xr -import cftime -import os -import collections +import os, types, json, cftime, collections ### single glacier output parent class ### @dataclass @@ -31,25 +29,21 @@ class single_glacier: """ glacier_rgi_table : pd.DataFrame dates_table : pd.DataFrame - wateryear : bool pygem_version : float - user_info : dict - outdir : str gcm_name : str scenario : str realization : str - calib_opt : int sim_iters : int modelprms : dict - ba_opt : int - gcm_bc_startyr : int - gcm_startyr : int - gcm_endyr: int + 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() @@ -64,32 +58,32 @@ def set_fn(self, fn=None): self.outfn += f'{self.scenario}_' if self.realization: self.outfn += f'{self.realization}_' - if self.calib_opt: - self.outfn += f'{self.calib_opt}_' + 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"]}' + 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{self.ba_opt}_' - yr0 = self.gcm_bc_startyr + self.outfn += f'ba{pygem_prms.option_bias_adjustment}_' + yr0 = self.gcm_bc_startyear else: self.outfn += 'ba0_' - yr0 = self.gcm_startyr - if self.calib_opt: + yr0 = self.gcm_startyear + if pygem_prms.option_calibration: self.outfn += 'SETS_' self.outfn += f'{yr0}_' - self.outfn += f'{self.gcm_endyr}_' + self.outfn += f'{self.gcm_endyear}_' def get_fn(self): return self.outfn def set_time_vals(self): - if self.wateryear == 'hydro': + 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 self.wateryear == 'calendar': + 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 self.wateryear == 'custom': + 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] @@ -98,7 +92,22 @@ def set_time_vals(self): self.year_values = np.concatenate((self.year_values, np.array([self.annual_columns[-1] + 1]))) def model_params_record(self): - pass + 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 @@ -179,9 +188,10 @@ def create_xr_ds(self): 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': self.user_info['institution'], - 'history': f'Created by {self.user_info["name"]} ({self.user_info["email"]}) on ' + datetime.today().strftime('%Y-%m-%d'), - 'references': 'doi:10.3389/feart.2019.00331 and doi:10.1017/jog.2019.91'} + '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 @@ -198,7 +208,6 @@ class glacierwide_stats(single_glacier): """ Single glacier-wide statistics dataset """ - extra_vars : bool 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) @@ -302,7 +311,7 @@ def update_dicts(self): 'temporal_resolution': 'monthly', 'comment': 'off-glacier runoff from area where glacier no longer exists'} - if self.extra_vars: + 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'] = { From dca52be11a0fd7119f81c91a2b4c1035e4cd5398 Mon Sep 17 00:00:00 2001 From: "Brandon S. Tober" Date: Fri, 26 Apr 2024 11:54:08 -0400 Subject: [PATCH 16/45] Update pygem install documentation for development mode (#41) * output cleanup, minor bug fixes * dev install documentation update --- docs/install_pygem.md | 13 +++++++------ pyproject.toml | 1 + setup.py | 21 ++++----------------- 3 files changed, 12 insertions(+), 23 deletions(-) mode change 100755 => 100644 setup.py 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/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 From 81dd27afb7721b863d919eb8a4d4b7dba7349905 Mon Sep 17 00:00:00 2001 From: "Brandon S. Tober" Date: Wed, 3 Apr 2024 10:33:37 -0400 Subject: [PATCH 17/45] Output class objects (#39) * bug fix in pygem_modelsetup.datesmodelrun(). pandas.date_range() has a timestamp range limitation based on the chosen resolution. The deafult is nanoseconds, which was reduced to seconds: https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#timeseries-timestamp-limits * check for differing climate data file paths or names based * framework for PyGEM output classes * added conda env file * output dataclasses * output cleanup * output binned stats * model setup group by thousands * output classes progress * Delete environment.yml --------- Co-authored-by: David Rounce --- pygem/output.py | 651 ++++++++++++++++++++++++++++++++++++++ pygem/pygem_modelsetup.py | 13 +- 2 files changed, 659 insertions(+), 5 deletions(-) create mode 100644 pygem/output.py diff --git a/pygem/output.py b/pygem/output.py new file mode 100644 index 00000000..1decf628 --- /dev/null +++ b/pygem/output.py @@ -0,0 +1,651 @@ +#!/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 cftime +import os +import 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 + wateryear : bool + pygem_version : float + user_info : dict + outdir : str + gcm_name : str + scenario : str + realization : str + calib_opt : int + sim_iters : int + modelprms : dict + ba_opt : int + gcm_bc_startyr : int + gcm_startyr : int + gcm_endyr: 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.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 self.calib_opt: + self.outfn += f'{self.calib_opt}_' + 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{self.ba_opt}_' + yr0 = self.gcm_bc_startyr + else: + self.outfn += 'ba0_' + yr0 = self.gcm_startyr + if self.calib_opt: + self.outfn += 'SETS_' + self.outfn += f'{yr0}_' + self.outfn += f'{self.gcm_endyr}_' + + def get_fn(self): + return self.outfn + + def set_time_vals(self): + if self.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 self.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 self.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): + pass + + 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': self.user_info['institution'], + 'history': f'Created by {self.user_info["name"]} ({self.user_info["email"]}) on ' + datetime.today().strftime('%Y-%m-%d'), + 'references': 'doi:10.3389/feart.2019.00331 and doi:10.1017/jog.2019.91'} + + 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 + """ + extra_vars : bool + + 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 self.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 50ecff90..ff6230bf 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -264,7 +264,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. @@ -319,7 +320,8 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' # 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]) @@ -327,8 +329,9 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' 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] @@ -502,7 +505,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() From 237729aa3f8092f12f2cbc6e526dddf47b34a9fc Mon Sep 17 00:00:00 2001 From: "Brandon S. Tober" Date: Wed, 3 Apr 2024 20:40:54 -0400 Subject: [PATCH 18/45] output cleanup, minor bug fixes (#40) --- pygem/output.py | 65 ++++++++++++++++++++++++++++--------------------- 1 file changed, 37 insertions(+), 28 deletions(-) diff --git a/pygem/output.py b/pygem/output.py index 1decf628..c499957f 100644 --- a/pygem/output.py +++ b/pygem/output.py @@ -19,9 +19,7 @@ import numpy as np import pandas as pd import xarray as xr -import cftime -import os -import collections +import os, types, json, cftime, collections ### single glacier output parent class ### @dataclass @@ -31,25 +29,21 @@ class single_glacier: """ glacier_rgi_table : pd.DataFrame dates_table : pd.DataFrame - wateryear : bool pygem_version : float - user_info : dict - outdir : str gcm_name : str scenario : str realization : str - calib_opt : int sim_iters : int modelprms : dict - ba_opt : int - gcm_bc_startyr : int - gcm_startyr : int - gcm_endyr: int + 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() @@ -64,32 +58,32 @@ def set_fn(self, fn=None): self.outfn += f'{self.scenario}_' if self.realization: self.outfn += f'{self.realization}_' - if self.calib_opt: - self.outfn += f'{self.calib_opt}_' + 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"]}' + 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{self.ba_opt}_' - yr0 = self.gcm_bc_startyr + self.outfn += f'ba{pygem_prms.option_bias_adjustment}_' + yr0 = self.gcm_bc_startyear else: self.outfn += 'ba0_' - yr0 = self.gcm_startyr - if self.calib_opt: + yr0 = self.gcm_startyear + if pygem_prms.option_calibration: self.outfn += 'SETS_' self.outfn += f'{yr0}_' - self.outfn += f'{self.gcm_endyr}_' + self.outfn += f'{self.gcm_endyear}_' def get_fn(self): return self.outfn def set_time_vals(self): - if self.wateryear == 'hydro': + 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 self.wateryear == 'calendar': + 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 self.wateryear == 'custom': + 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] @@ -98,7 +92,22 @@ def set_time_vals(self): self.year_values = np.concatenate((self.year_values, np.array([self.annual_columns[-1] + 1]))) def model_params_record(self): - pass + 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 @@ -179,9 +188,10 @@ def create_xr_ds(self): 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': self.user_info['institution'], - 'history': f'Created by {self.user_info["name"]} ({self.user_info["email"]}) on ' + datetime.today().strftime('%Y-%m-%d'), - 'references': 'doi:10.3389/feart.2019.00331 and doi:10.1017/jog.2019.91'} + '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 @@ -198,7 +208,6 @@ class glacierwide_stats(single_glacier): """ Single glacier-wide statistics dataset """ - extra_vars : bool 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) @@ -302,7 +311,7 @@ def update_dicts(self): 'temporal_resolution': 'monthly', 'comment': 'off-glacier runoff from area where glacier no longer exists'} - if self.extra_vars: + 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'] = { From a01cafbab24709bc69b17a64ec11cd6e6b0b87d7 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Thu, 7 Dec 2023 09:30:31 +0100 Subject: [PATCH 19/45] scale the calving-k(yield strength) --- pygem/glacierdynamics.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index c79132e1..8990f2d6 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -167,7 +167,11 @@ def balance_thickness(yield_strength, bed_elev): else: D = 0 return (2 * yield_strength / (RHO_ICE * G)) + np.sqrt( +<<<<<<< HEAD (RHO_SEA * (D ** 2) / RHO_ICE) + ((2 * yield_strength / (RHO_ICE * G)) ** 2)) +======= + (RHO_SEA * (D ** 2) / RHO_ICE) + ((2 * yield_strength / (RHO_ICE * G)) ** 2)) +>>>>>>> cc18294 (scale the calving-k(yield strength)) # TODO: Check on exponent on last term. In Ultee & Bassis 2016, this is squared, but in Ultee & Bassis 2020 supplement, it isn't. # --------------------------------------------------------------------------- From 367a9b906a9b64e529343655e3f14ba259800b2f Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Wed, 8 May 2024 15:41:42 +0200 Subject: [PATCH 20/45] add the try-except part in the function get_annual_mb, to export the traceback error --- pygem/massbalance.py | 965 ++++++++++++++++++++++--------------------- 1 file changed, 496 insertions(+), 469 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 190c7055..a30a48e2 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 @@ -49,7 +50,7 @@ def __init__(self, gdir, modelprms, glacier_rgi_table, switch to run the model in reverse or not (may be irrelevant after converting to OGGM's setup) """ - print("PYGEM: flowlines") + print("the input flowlines in PyGEMMassBalance Model is: ") print(type(fls).__name__) if debug: @@ -197,7 +198,7 @@ def get_monthly_mb(self, heights, year=None, fls=None, fl_id=None, return mb - def get_annual_mb(self, heights, year=None, fls=None, fl_id=None, + def get_annual_mb(self, heights, year=None, fls = None, fl_id = 0, debug=False, option_areaconstant=False, year_month=None): """FIXED FORMAT FOR THE FLOWLINE MODEL @@ -215,486 +216,512 @@ 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 - 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 + #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()) + 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] - 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] - 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: - 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] + 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] - 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: - seconds_in_month = self.dayspermonth[12*year+round(year_month*12)]* 24 * 3600 - mb = (self.glac_bin_massbalclim[:,12*year+round(year_month*12)] - * pygem_prms.density_water / pygem_prms.density_ice - /seconds_in_month) - print("index for mb") - print(12*year+int(year_month*12)) - print(year_month) - print(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+round(year_month*12)]* 24 * 3600 + mb = (self.glac_bin_massbalclim[:,12*year+round(year_month*12)] + * pygem_prms.density_water / pygem_prms.density_ice + /seconds_in_month) + print("index for mb") + print(12*year+int(year_month*12)) + print(year_month) + print(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 + except: + print(traceback.format_exc()) #%% def _convert_glacwide_results(self, year, glacier_area, heights, From 8984c2ea6b902af58d323ef63e9013fa8a2e9105 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Wed, 8 May 2024 20:51:52 +0200 Subject: [PATCH 21/45] modify the equation in the function balance_thickness from (RHO_SEA * (D ** 2) / RHO_ICE) + ((2 * yield_strength / (RHO_ICE * G)) ** 1)) to (RHO_SEA * (D ** 2) / RHO_ICE) + ((2 * yield_strength / (RHO_ICE * G)) ** 2)) --- pygem/glacierdynamics.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index 8990f2d6..c79132e1 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -167,11 +167,7 @@ def balance_thickness(yield_strength, bed_elev): else: D = 0 return (2 * yield_strength / (RHO_ICE * G)) + np.sqrt( -<<<<<<< HEAD (RHO_SEA * (D ** 2) / RHO_ICE) + ((2 * yield_strength / (RHO_ICE * G)) ** 2)) -======= - (RHO_SEA * (D ** 2) / RHO_ICE) + ((2 * yield_strength / (RHO_ICE * G)) ** 2)) ->>>>>>> cc18294 (scale the calving-k(yield strength)) # TODO: Check on exponent on last term. In Ultee & Bassis 2016, this is squared, but in Ultee & Bassis 2020 supplement, it isn't. # --------------------------------------------------------------------------- From eb1bcebc40eb5d91944f7dea69b287d88b0479a0 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Tue, 25 Jun 2024 20:51:18 +0200 Subject: [PATCH 22/45] revise the seconds_in_month = self.dayspermonth[12*year+int(year_month*12)]* 24 * 3600 to make sure the float year is transfered to the correct int year and month, which was self.dayspermonth[12*year+round(year_month*12)]* 24 * 3600, in which the last one will be out of the range, same to mb, in the function get_annual_mb; changed the default fl_id = -1 from 0, actually it doesn't matter, because the elevation band flowline only include one flowline, but the -1 always refers to the main flowline(longest) if they includes multiple flowlines; comment some print hints --- pygem/massbalance.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index a30a48e2..1edaec74 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -198,7 +198,7 @@ def get_monthly_mb(self, heights, year=None, fls=None, fl_id=None, return mb - def get_annual_mb(self, heights, year=None, fls = None, fl_id = 0, + 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 @@ -679,8 +679,8 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = 0, 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+round(year_month*12)]* 24 * 3600 - mb = (self.glac_bin_massbalclim[:,12*year+round(year_month*12)] + seconds_in_month = self.dayspermonth[12*year+int(year_month*12)]* 24 * 3600 + 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") From 7271dc618351d4928b1add58bb6d1ebbdcd2627a Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Tue, 25 Jun 2024 20:53:08 +0200 Subject: [PATCH 23/45] change the parameter cfg.PARAMS['continue_on_error'] = False, for the single glacier test --- pygem/oggm_compat.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygem/oggm_compat.py b/pygem/oggm_compat.py index 62e4a3d8..124ccd3d 100755 --- a/pygem/oggm_compat.py +++ b/pygem/oggm_compat.py @@ -50,7 +50,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 @@ -143,7 +143,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 From 345924b282afb3b2ce5a7cd01c6f9107be2c5bb3 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Tue, 25 Jun 2024 20:54:54 +0200 Subject: [PATCH 24/45] import the traceback to trace the error and add some print hints --- pygem/pygem_modelsetup.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pygem/pygem_modelsetup.py b/pygem/pygem_modelsetup.py index ff6230bf..67fab4c6 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 @@ -306,7 +307,7 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' 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): @@ -317,6 +318,7 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' csv_regionO1 = pd.read_csv(rgi_fp + rgi_fn) except: csv_regionO1 = pd.read_csv(rgi_fp + rgi_fn, encoding='latin1') + print(traceback.format_exc()) # Populate glacer_table with the glaciers of interest if rgi_regionsO2 == 'all' and rgi_glac_number == 'all': @@ -354,7 +356,7 @@ 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.) From 37bb1841ae6a474932d2e2bf906f61ed7beec9fd Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Tue, 25 Jun 2024 21:00:25 +0200 Subject: [PATCH 25/45] comment some print hints --- pygem/massbalance.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 1edaec74..5e7dac9b 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -50,8 +50,8 @@ def __init__(self, gdir, modelprms, glacier_rgi_table, 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__) + #print("the input flowlines in PyGEMMassBalance Model is: ") + #print(type(fls).__name__) if debug: print('\n\nDEBUGGING MASS BALANCE FUNCTION\n\n') @@ -218,14 +218,14 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, """ #year = int(year) year=int(year) - print("Year is: ", year, "Month is:", year_month) + #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("#########################") + # 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") @@ -237,7 +237,7 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, # print(traceback.format_exc()) if (fls is not None) and (fl_id is not None): - print("fls and fl_id are not None") + #print("fls and fl_id are not None") try: fl = fls[fl_id] except: @@ -271,7 +271,7 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, # Glacier indices glac_idx_t0 = glacier_area_t0.nonzero()[0] - print('the indices of glacier nonzero at t0 is:',glac_idx_t0) + #print('the indices of glacier nonzero at t0 is:',glac_idx_t0) nbins = heights.shape[0] nmonths = self.glacier_gcm_temp.shape[0] From 4bd65412b2c5158dcb3d308590a86a13566b2774 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 22 Jul 2024 15:08:17 +0200 Subject: [PATCH 26/45] add a try/except part in the def of single_flowline_glacier_directory --- pygem/oggm_compat.py | 49 ++++++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/pygem/oggm_compat.py b/pygem/oggm_compat.py index 124ccd3d..8a1a3ea9 100755 --- a/pygem/oggm_compat.py +++ b/pygem/oggm_compat.py @@ -3,6 +3,7 @@ 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 @@ -79,29 +80,33 @@ 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] - - # 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) + try: + # 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] - gdir = gdirs[0] + # 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] + except: + print("Somenthing wrong with the single_flowline_glacier_directory") + print(traceback.format_exc()) return gdir From cbdaf93954def7a9a0dd4fa591e6d7ccef50bc3e Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 19 Aug 2024 10:17:34 +0200 Subject: [PATCH 27/45] revising the print hints --- pygem/massbalance.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 5e7dac9b..51b0dece 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -684,9 +684,9 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, * pygem_prms.density_water / pygem_prms.density_ice /seconds_in_month) print("index for mb") - print(12*year+int(year_month*12)) - print(year_month) - print(year) + 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) From 0fc5010780ca8d32b5876b1030cae855269872a1 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Tue, 27 Aug 2024 15:02:59 +0200 Subject: [PATCH 28/45] add a print hint about of the off glacier idx --- pygem/massbalance.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 51b0dece..7fac51c9 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -298,7 +298,8 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, 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] - print('****************** off glacier idx is*******************',offglac_idx) + #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 From a8ef5bff5510bbe6236e189682195ea2005d14d3 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Tue, 3 Sep 2024 17:56:31 +0200 Subject: [PATCH 29/45] add print hint of year and month in def get_monthly_mb --- pygem/massbalance.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 7fac51c9..aa5fa92c 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -192,7 +192,8 @@ def get_monthly_mb(self, heights, year=None, fls=None, fl_id=None, 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) From 80571ec215491276f97db74ee2a95682290d1132 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Wed, 18 Sep 2024 14:38:49 +0200 Subject: [PATCH 30/45] add print hints --- pygem/massbalance.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index aa5fa92c..2b6fc595 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -680,8 +680,9 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -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) + 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) @@ -720,7 +721,7 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, # plt.show() # # print('mb_filled:', mb_filled) - + print("**************** get_annual_mb end ****************") return mb_filled except: print(traceback.format_exc()) From 5e2685260c19c19e79e6566b91eaac500a7fbe66 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Thu, 24 Oct 2024 14:19:47 +0200 Subject: [PATCH 31/45] add the glac_length_change as extra variable to store; comment one print hint --- pygem/massbalance.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 2b6fc595..0f3fea87 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -144,7 +144,7 @@ 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_wide_length_annual = np.zeros(self.nyears+1) + self.glac_length_change = np.zeros(self.nmonths) self.glac_wide_massbaltotal = np.zeros(self.nmonths) self.glac_wide_runoff = np.zeros(self.nmonths) self.glac_wide_snowline = np.zeros(self.nmonths) @@ -276,7 +276,7 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, nbins = heights.shape[0] nmonths = self.glacier_gcm_temp.shape[0] - print('nbins is',nbins,'nmonths is :',nmonths) + #print('nbins is',nbins,'nmonths is :',nmonths) # Local variables bin_precsnow = np.zeros((nbins,nmonths)) From 66b1fab3d1145232ff4f670bd858e5024a0bafc7 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Thu, 24 Oct 2024 14:26:48 +0200 Subject: [PATCH 32/45] add the millan22 thickness in the gdir --- pygem/oggm_compat.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/pygem/oggm_compat.py b/pygem/oggm_compat.py index 8a1a3ea9..ae6f8b63 100755 --- a/pygem/oggm_compat.py +++ b/pygem/oggm_compat.py @@ -13,6 +13,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): @@ -87,6 +90,9 @@ def single_flowline_glacier_directory(rgi_id, reset=pygem_prms.overwrite_gdirs, 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') + # add the millan22 thickness and velocity + #workflow.execute_entity_task(millan22.velocity_to_gdir, gdirs) + workflow.execute_entity_task(millan22.thickness_to_gdir, gdirs) # Compute all the stuff list_tasks = [ @@ -186,6 +192,9 @@ def single_flowline_glacier_directory_with_calving(rgi_id, reset=pygem_prms.over if not gdirs[0].is_tidewater: raise ValueError(f'{rgi_id} is not tidewater!') + # add the millan22 thickness and velocity + #workflow.execute_entity_task(millan22.velocity_to_gdir, gdirs) + workflow.execute_entity_task(millan22.thickness_to_gdir, gdirs) # Compute all the stuff list_tasks = [ # Consensus ice thickness From 54d19ffc5e6ca968548103c3ee71d4e995e200ef Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 4 Nov 2024 16:40:42 +0100 Subject: [PATCH 33/45] in the def of ensure_mass_conservation, add the new argument Dynamic_step_Monthly,to control the annual/monthly dynamic run --- pygem/massbalance.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 0f3fea87..2c849a3d 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -876,7 +876,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 @@ -889,14 +889,18 @@ 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) #%% Dynamic running step - # if the dynamic step is monthly, the volume is monthly, need to be calculated as annual - Dynamic_step_Monthly =True + #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) From a1aebfe78f81d3d40127e2b0d644908f4520516c Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Wed, 27 Nov 2024 21:05:52 +0100 Subject: [PATCH 34/45] add a input argument k_calving_str in the function single_flowline_glacier_directory_with_calving, to generate the path of oggm_gdir for the paralle computing for assemble parts related with data assimilation --- pygem/oggm_compat.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygem/oggm_compat.py b/pygem/oggm_compat.py index ae6f8b63..fdca62a1 100755 --- a/pygem/oggm_compat.py +++ b/pygem/oggm_compat.py @@ -119,7 +119,7 @@ def single_flowline_glacier_directory(rgi_id, reset=pygem_prms.overwrite_gdirs, 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) @@ -167,7 +167,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: From b7d840f76c0e669736aa7d504ea8f42d4f8f68be Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 10 Mar 2025 14:20:20 +0100 Subject: [PATCH 35/45] add a print hint of mb_clim_cn --- pygem/shop/mbdata.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pygem/shop/mbdata.py b/pygem/shop/mbdata.py index 84b51a10..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] From 8a9a2adb95ce91eea303ab79753acd268c1589f6 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 10 Mar 2025 14:22:07 +0100 Subject: [PATCH 36/45] add the glac_wide_massbalclim as avariable in class PyGEMMassBalance --- pygem/massbalance.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 2c849a3d..55cec97f 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -146,6 +146,7 @@ def __init__(self, gdir, modelprms, glacier_rgi_table, 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) @@ -798,6 +799,9 @@ 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_massbaltotal[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: From 9bff23c9e2fbcd62763e273dc57801ab682b280e Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 10 Mar 2025 14:29:20 +0100 Subject: [PATCH 37/45] modify the functions sinlge_flowline_glacier_directory and single_flowline_glacier_directory_with_calving, the contents are the same just with better workflow --- pygem/oggm_compat.py | 66 +++++++++++++++++++++++--------------------- 1 file changed, 34 insertions(+), 32 deletions(-) diff --git a/pygem/oggm_compat.py b/pygem/oggm_compat.py index fdca62a1..eb395d43 100755 --- a/pygem/oggm_compat.py +++ b/pygem/oggm_compat.py @@ -1,5 +1,6 @@ """ PYGEM-OGGGM COMPATIBILITY FUNCTIONS """ # External libraries +import os import numpy as np import pandas as pd import netCDF4 @@ -73,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 @@ -88,28 +91,28 @@ def single_flowline_glacier_directory(rgi_id, reset=pygem_prms.overwrite_gdirs, 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') + 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) - - # Compute all the stuff - list_tasks = [ - # Consensus ice thickness - icethickness.consensus_gridded, - # Mass balance data - mbdata.mb_df_to_gdir] - + #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 + # debris thickness and melt enhancement factors 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] + 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()) @@ -186,27 +189,26 @@ 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] - if not gdirs[0].is_tidewater: + if not gdir.is_tidewater: raise ValueError(f'{rgi_id} is not tidewater!') # add the millan22 thickness and velocity #workflow.execute_entity_task(millan22.velocity_to_gdir, gdirs) - workflow.execute_entity_task(millan22.thickness_to_gdir, gdirs) - # Compute all the stuff - list_tasks = [ - # Consensus ice thickness - icethickness.consensus_gridded, - # Mass balance data - mbdata.mb_df_to_gdir] + #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): From 3b38d7fabd6eba8d6af079329bba50badfea23d1 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Tue, 25 Mar 2025 10:25:20 +0100 Subject: [PATCH 38/45] add the input argument filesuffix for the function debris_binned, and adding it when calling gdir.read_pickle,and gdir.write_pickle --- pygem/shop/debris.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) 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 From a12936874959e514e9b8132c98c61a38cfe505e5 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Wed, 2 Apr 2025 16:09:18 +0200 Subject: [PATCH 39/45] changed the print hints for the csv file encoding with 'latin1' in the def of selectglaciersrgitable --- pygem/pygem_modelsetup.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygem/pygem_modelsetup.py b/pygem/pygem_modelsetup.py index 67fab4c6..8f24dce1 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -318,7 +318,8 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' csv_regionO1 = pd.read_csv(rgi_fp + rgi_fn) except: csv_regionO1 = pd.read_csv(rgi_fp + rgi_fn, encoding='latin1') - print(traceback.format_exc()) + print('..latin1 encoding..') + #print(traceback.format_exc()) # Populate glacer_table with the glaciers of interest if rgi_regionsO2 == 'all' and rgi_glac_number == 'all': From ad3976b823152a01c2a209c2a3951631e99430c0 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Wed, 16 Apr 2025 14:36:53 +0200 Subject: [PATCH 40/45] add the print hints for the assertion of the height = fl.surfac_h in the function get_annual_mb --- pygem/massbalance.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 55cec97f..bfc8f11a 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -247,6 +247,13 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, 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 @@ -289,7 +296,7 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, refreeze_potential = np.zeros(nbins) if self.glacier_area_initial.sum() > 0: - # if len(glac_idx_t0) > 0: + # if len(glac_idx_t0) > 0: # Surface type [0=off-glacier, 1=ice, 2=snow, 3=firn, 4=debris] if year == 0: @@ -664,15 +671,15 @@ def get_annual_mb(self, heights, year=None, fls = None, fl_id = -1, 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]) + ## 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: From 2aaf7af54dc256f989b9105bf9af9dec0c901fbe Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Thu, 17 Apr 2025 15:05:50 +0200 Subject: [PATCH 41/45] revised the self.glac_wide_massbalclim as acc+refreeze-melt, under the def of _convert_glacwide_results --- pygem/massbalance.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index bfc8f11a..4badbbb5 100755 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -808,7 +808,8 @@ def _convert_glacwide_results(self, year, glacier_area, heights, - 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_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)]) # 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: From c4b4e9f1f08e740c0933710774bde642d98a5fb9 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Mon, 21 Apr 2025 21:11:14 +0200 Subject: [PATCH 42/45] add the hardcode gdir.is_tidewater = True, in the single_flowline_glacier_directory_with_calving, because the RGI DATASET DOESNOT HAVE THE CORRECT TERMINUS TYPE, WE ASSUME GALCIERS WITH FRONATAL ABLATION DATASET ARE ALL TIDEWATER GLACIER --- pygem/oggm_compat.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pygem/oggm_compat.py b/pygem/oggm_compat.py index eb395d43..7cd94367 100755 --- a/pygem/oggm_compat.py +++ b/pygem/oggm_compat.py @@ -192,6 +192,9 @@ def single_flowline_glacier_directory_with_calving(rgi_id, reset=pygem_prms.over 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 gdir.is_tidewater: raise ValueError(f'{rgi_id} is not tidewater!') From e31145cdbc0f47ef515a69e59686d0bc9932b248 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Wed, 6 Aug 2025 19:19:39 +0200 Subject: [PATCH 43/45] comment two print hints --- pygem/pygem_modelsetup.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygem/pygem_modelsetup.py b/pygem/pygem_modelsetup.py index 8f24dce1..089a6b50 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -300,7 +300,7 @@ 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('rfp',os.listdir(rgi_fp)) print('rgi_glc_num',rgi_glac_number) glacier_table = pd.DataFrame() for region in rgi_regionsO1: @@ -313,7 +313,7 @@ def selectglaciersrgitable(glac_no=None, rgi_regionsO1=None, rgi_regionsO2='all' for i in os.listdir(rgi_fp): if i.startswith(str(region).zfill(2)) and i.endswith('.csv'): rgi_fn = i - print('regs:',i) + #print('regs:',i) try: csv_regionO1 = pd.read_csv(rgi_fp + rgi_fn) except: From 8b95a2a0e37d2ca1d5487b9fc9653d9796ca32a4 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Thu, 14 Aug 2025 10:04:54 +0200 Subject: [PATCH 44/45] add detailed value info in the assert error message of gcm_prec_biasadj_fra in the function prec_biasadj_opt1 --- pygem/gcmbiasadj.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygem/gcmbiasadj.py b/pygem/gcmbiasadj.py index 6419e54e..576c1674 100755 --- a/pygem/gcmbiasadj.py +++ b/pygem/gcmbiasadj.py @@ -322,7 +322,8 @@ def prec_biasadj_opt1(ref_prec, ref_elev, gcm_prec, dates_table_ref, dates_table 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') + 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' From 02e8b5c58e1e540102d2390b308941f390fa5de0 Mon Sep 17 00:00:00 2001 From: Ruitang Yang Date: Thu, 14 Aug 2025 12:18:48 +0200 Subject: [PATCH 45/45] adjust the threshold of gcm_prec_biasadj_frac as 0.45 from 0.5 (lower) --- pygem/gcmbiasadj.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygem/gcmbiasadj.py b/pygem/gcmbiasadj.py index 576c1674..2a2d6892 100755 --- a/pygem/gcmbiasadj.py +++ b/pygem/gcmbiasadj.py @@ -321,7 +321,7 @@ 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, ( + 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'