diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite.yml index 4bb4a542..2c30c1bb 100644 --- a/.github/workflows/test_suite.yml +++ b/.github/workflows/test_suite.yml @@ -73,4 +73,5 @@ jobs: python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_01_basics.py python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_02_config.py python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_03_notebooks.py - python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_04_postproc.py + python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_04_auxiliary.py + python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_05_postproc.py diff --git a/pygem/bin/postproc/postproc_binned_subannual_thick.py b/pygem/bin/postproc/postproc_binned_subannual_thick.py index 0af94963..04e0c93f 100644 --- a/pygem/bin/postproc/postproc_binned_subannual_thick.py +++ b/pygem/bin/postproc/postproc_binned_subannual_thick.py @@ -333,12 +333,10 @@ def main(): if args.simpath: # filter out non-file paths simpath = [p for p in args.simpath if os.path.isfile(p)] - elif args.simdir: # get list of sims simpath = sorted(glob.glob(args.simdir + '/*.nc')) if simpath: - print(simpath) # number of cores for parallel processing if args.ncores > 1: ncores = int(np.min([len(simpath), args.ncores])) diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py index 069abc93..43fb2285 100755 --- a/pygem/bin/run/run_simulation.py +++ b/pygem/bin/run/run_simulation.py @@ -986,7 +986,9 @@ def run(list_packed_vars): if debug: fig, ax = plt.subplots(1) - graphics.plot_modeloutput_section(ev_model, ax=ax) + graphics.plot_modeloutput_section( + ev_model, ax=ax, lnlabel=f'Glacier year {args.sim_startyear}' + ) diag = ev_model.run_until_and_store(args.sim_endyear + 1) ev_model.mb_model.glac_wide_volume_annual[-1] = diag.volume_m3[-1] @@ -1056,9 +1058,10 @@ def run(list_packed_vars): ) if debug: - print('New glacier vol', ev_model.volume_m3) - graphics.plot_modeloutput_section(ev_model) - plt.show() + fig, ax = plt.subplots(1) + graphics.plot_modeloutput_section( + ev_model, ax=ax, lnlabel=f'Glacier year {args.sim_startyear}' + ) _, diag = ev_model.run_until_and_store(args.sim_endyear + 1) # print('shape of volume:', ev_model.mb_model.glac_wide_volume_annual.shape, diag.volume_m3.shape) @@ -1113,25 +1116,20 @@ def run(list_packed_vars): years = np.arange(args.sim_startyear, args.sim_endyear + 1) mb_all = [] for year in years: - mb_annual = mbmod.get_annual_mb( + # Calculate annual mass balance + mbmod.get_annual_mb( nfls[0].surface_h, fls=nfls, fl_id=0, year=year, - debug=True, + debug=debug, ) - mb_mwea = ( - mb_annual - * 365 - * 24 - * 3600 - * pygem_prms['constants']['density_ice'] - / pygem_prms['constants']['density_water'] + # Record glacierwide annual mass balance + t_start, t_stop = mbmod.get_step_inds(year) + mb_all.append( + mbmod.glac_wide_massbaltotal[t_start : t_stop + 1].sum() + / mbmod.glacier_area_initial.sum() ) - glac_wide_mb_mwea = ( - mb_mwea * mbmod.glacier_area_initial - ).sum() / mbmod.glacier_area_initial.sum() - mb_all.append(glac_wide_mb_mwea) mbmod.glac_wide_area_annual[-1] = mbmod.glac_wide_area_annual[0] mbmod.glac_wide_volume_annual[-1] = mbmod.glac_wide_volume_annual[0] diag['area_m2'] = mbmod.glac_wide_area_annual @@ -1149,15 +1147,13 @@ def run(list_packed_vars): np.round(np.median(mb_all), 3), ) - # mb_em_mwea = run_emulator_mb(modelprms) - # print(' emulator mb:', np.round(mb_em_mwea,3)) - # mb_em_sims.append(mb_em_mwea) - # Record output for successful runs if successful_run: if args.option_dynamics is not None: if debug: - graphics.plot_modeloutput_section(ev_model, ax=ax, srfls='--') + graphics.plot_modeloutput_section( + ev_model, ax=ax, srfls='--', lnlabel=f'Glacier year {args.sim_endyear + 1}' + ) plt.figure() diag.volume_m3.plot() plt.show() @@ -1192,7 +1188,7 @@ def run(list_packed_vars): print(' mb_mbmod [mwea]:', np.round(mb_mwea_mbmod, 2)) if np.abs(mb_mwea_diag - mb_mwea_mbmod) > 1e-6: - ev_model.mb_model.ensure_mass_conservation(diag, dates_table) + ev_model.mb_model.ensure_mass_conservation(diag) if debug: print( diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index b4ca6bdb..16f2e910 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -48,7 +48,6 @@ def __init__( inplace=False, debug=True, option_areaconstant=False, - spinupyears=0, constantarea_years=0, **kwargs, ): @@ -87,12 +86,10 @@ def __init__( ) self.option_areaconstant = option_areaconstant self.constantarea_years = constantarea_years - self.spinupyears = spinupyears self.glac_idx_initial = [fl.thick.nonzero()[0] for fl in flowlines] - self.y0 = 0 + self.y0 = y0 self.is_tidewater = is_tidewater self.water_level = water_level - # widths_t0 = flowlines[0].widths_m # area_v1 = widths_t0 * flowlines[0].dx_meter # print('area v1:', area_v1.sum()) @@ -326,8 +323,13 @@ def run_until_and_store(self, y1, run_path=None, diag_path=None, store_monthly_s def updategeometry(self, year, debug=False): """Update geometry for a given year""" + # get year index + year_idx = self.mb_model.get_year_index(year) + # get time step indices - note indexing should be [t_start:t_stop+1] to include final step in year + t_start, t_stop = self.mb_model.get_step_inds(year) if debug: - print('year:', year) + print('year:', year, f'({year_idx})') + print('time steps:', f'[{t_start}, {t_stop}]') # Loop over flowlines for fl_id, fl in enumerate(self.fls): @@ -338,8 +340,8 @@ def updategeometry(self, year, debug=False): width_t0 = self.fls[fl_id].widths_m.copy() # CONSTANT AREAS - # Mass redistribution ignored for calibration and spinup years (glacier properties constant) - if (self.option_areaconstant) or (year < self.spinupyears) or (year < self.constantarea_years): + # Mass redistribution ignored for calibration years (glacier properties constant) + if (self.option_areaconstant) or (year < self.y0 + self.constantarea_years): # run mass balance glac_bin_massbalclim_annual = self.mb_model.get_annual_mb( heights, fls=self.fls, fl_id=fl_id, year=year, debug=False @@ -380,7 +382,7 @@ def updategeometry(self, year, debug=False): # If frontal ablation more than bin volume, remove entire bin if fa_m3 > vol_last: # Record frontal ablation (m3 w.e.) in mass balance model for output - self.mb_model.glac_bin_frontalablation[last_idx, int(12 * (year + 1) - 1)] = ( + self.mb_model.glac_bin_frontalablation[last_idx, t_stop + 1] = ( vol_last * pygem_prms['constants']['density_ice'] / pygem_prms['constants']['density_water'] @@ -397,7 +399,7 @@ def updategeometry(self, year, debug=False): section_t0[last_idx] = section_t0[last_idx] - fa_m3 / fl.dx_meter self.fls[fl_id].section = section_t0 # Record frontal ablation(m3 w.e.) - self.mb_model.glac_bin_frontalablation[last_idx, int(12 * (year + 1) - 1)] = ( + self.mb_model.glac_bin_frontalablation[last_idx, t_stop + 1] = ( fa_m3 * pygem_prms['constants']['density_ice'] / pygem_prms['constants']['density_water'] @@ -433,11 +435,8 @@ def updategeometry(self, year, debug=False): heights, fls=self.fls, fl_id=fl_id, year=year, debug=False ) sec_in_year = ( - self.mb_model.dates_table.loc[12 * year : 12 * (year + 1) - 1, 'daysinmonth'].values.sum() - * 24 - * 3600 + self.mb_model.dates_table.iloc[t_start : t_stop + 1]['days_in_step'].values.sum() * 24 * 3600 ) - # print(' volume change [m3]:', (glac_bin_massbalclim_annual * sec_in_year * # (width_t0 * fl.dx_meter)).sum()) # print(glac_bin_masssbalclim_annual) @@ -469,14 +468,13 @@ def updategeometry(self, year, debug=False): # Record glacier properties (volume [m3], area [m2], thickness [m], width [km]) # record the next year's properties as well # 'year + 1' used so the glacier properties are consistent with mass balance computations - year = int(year) # required to ensure proper indexing with run_until_and_store (10/21/2020) glacier_area = fl.widths_m * fl.dx_meter glacier_area[fl.thick == 0] = 0 - self.mb_model.glac_bin_area_annual[:, year + 1] = glacier_area - self.mb_model.glac_bin_icethickness_annual[:, year + 1] = fl.thick - self.mb_model.glac_bin_width_annual[:, year + 1] = fl.widths_m - self.mb_model.glac_wide_area_annual[year + 1] = glacier_area.sum() - self.mb_model.glac_wide_volume_annual[year + 1] = (fl.section * fl.dx_meter).sum() + self.mb_model.glac_bin_area_annual[:, year_idx + 1] = glacier_area + self.mb_model.glac_bin_icethickness_annual[:, year_idx + 1] = fl.thick + self.mb_model.glac_bin_width_annual[:, year_idx + 1] = fl.widths_m + self.mb_model.glac_wide_area_annual[year_idx + 1] = glacier_area.sum() + self.mb_model.glac_wide_volume_annual[year_idx + 1] = (fl.section * fl.dx_meter).sum() # %% ----- FRONTAL ABLATION ----- def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, calving_k=None, debug=False): diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 51482ceb..fa4a010c 100644 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -210,6 +210,18 @@ def get_year_index(self, year): assert year in self.years, f'{year} not found in model dates table' return self.year_to_index[year] + def get_step_inds(self, year): + """ + Get time step start and stop indices for the specified model year. + + Both indices are inclusive, so standard Python slicing would use + [t_start : t_stop + 1]. + """ + step_idxs = np.where(self.dates_table.year == int(year))[0] + if step_idxs.size == 0: + raise ValueError(f'Year {year} not found in dates_table.') + return step_idxs[0], step_idxs[-1] + def get_annual_mb( self, heights, @@ -243,12 +255,8 @@ def get_annual_mb( # get year index year_idx = self.get_year_index(year) - step_idxs = np.where(self.dates_table.year == int(year)) - # get start step for 0th step of specified year - t_start = step_idxs[0][0] - # get stop step for specified year - # note, this is 1 greater than the final step which to include - python indexing will not include this step - t_stop = step_idxs[0][-1] + 1 + # get time step indices - note indexing should be [t_start:t_stop+1] to include final step in year + t_start, t_stop = self.get_step_inds(year) fl = fls[fl_id] @@ -303,11 +311,11 @@ def get_annual_mb( # AIR TEMPERATURE: Downscale the gcm temperature [deg C] to each bin # Downscale using gcm and glacier lapse rates # T_bin = T_gcm + lr_gcm * (z_ref - z_gcm) + lr_glac * (z_bin - z_ref) + tempchange - self.bin_temp[:, t_start:t_stop] = ( - self.glacier_gcm_temp[t_start:t_stop] - + self.glacier_gcm_lrgcm[t_start:t_stop] + self.bin_temp[:, t_start : t_stop + 1] = ( + self.glacier_gcm_temp[t_start : t_stop + 1] + + self.glacier_gcm_lrgcm[t_start : t_stop + 1] * (self.glacier_rgi_table.loc[pygem_prms['mb']['option_elev_ref_downscale']] - self.glacier_gcm_elev) - + self.glacier_gcm_lrglac[t_start:t_stop] + + self.glacier_gcm_lrglac[t_start : t_stop + 1] * (heights - self.glacier_rgi_table.loc[pygem_prms['mb']['option_elev_ref_downscale']])[:, np.newaxis] + self.modelprms['tbias'] ) @@ -315,8 +323,8 @@ def get_annual_mb( # PRECIPITATION/ACCUMULATION: Downscale the precipitation (liquid and solid) to each bin # Precipitation using precipitation factor and precipitation gradient # P_bin = P_gcm * prec_factor * (1 + prec_grad * (z_bin - z_ref)) - bin_precsnow[:, t_start:t_stop] = ( - self.glacier_gcm_prec[t_start:t_stop] + bin_precsnow[:, t_start : t_stop + 1] = ( + self.glacier_gcm_prec[t_start : t_stop + 1] * self.modelprms['kp'] * ( 1 @@ -345,8 +353,8 @@ def get_annual_mb( height_75 = heights[glac_idx_upper25].min() glac_idx_75 = np.where(heights == height_75)[0][0] # exponential decay - bin_precsnow[glac_idx_upper25, t_start:t_stop] = ( - bin_precsnow[glac_idx_75, t_start:t_stop] + bin_precsnow[glac_idx_upper25, t_start : t_stop + 1] = ( + bin_precsnow[glac_idx_75, t_start : t_stop + 1] * np.exp( -1 * (heights[glac_idx_upper25] - height_75) @@ -356,61 +364,63 @@ def get_annual_mb( # Precipitation cannot be less than 87.5% of the maximum accumulation elsewhere on the glacier # compute max values for each step over glac_idx_t0, compare, and replace if needed max_values = np.tile( - 0.875 * np.max(bin_precsnow[glac_idx_t0, t_start:t_stop], axis=0), + 0.875 * np.max(bin_precsnow[glac_idx_t0, t_start : t_stop + 1], axis=0), (8, 1), ) - uncorrected_values = bin_precsnow[glac_idx_upper25, t_start:t_stop] + uncorrected_values = bin_precsnow[glac_idx_upper25, t_start : t_stop + 1] corrected_values = np.max(np.stack([uncorrected_values, max_values], axis=0), axis=0) - bin_precsnow[glac_idx_upper25, t_start:t_stop] = corrected_values + bin_precsnow[glac_idx_upper25, t_start : t_stop + 1] = corrected_values # Separate total precipitation into liquid (bin_prec) and solid (bin_acc) if pygem_prms['mb']['option_accumulation'] == 1: # if temperature above threshold, then rain ( - self.bin_prec[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] > self.modelprms['tsnow_threshold'] + self.bin_prec[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] ] - ) = bin_precsnow[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] > self.modelprms['tsnow_threshold'] + ) = bin_precsnow[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] ] # if temperature below threshold, then snow ( - self.bin_acc[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] <= self.modelprms['tsnow_threshold'] + self.bin_acc[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] ] - ) = bin_precsnow[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] <= self.modelprms['tsnow_threshold'] + ) = bin_precsnow[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] ] elif pygem_prms['mb']['option_accumulation'] == 2: # if temperature between min/max, then mix of snow/rain using linear relationship between min/max - self.bin_prec[:, t_start:t_stop] = ( - 0.5 + (self.bin_temp[:, t_start:t_stop] - self.modelprms['tsnow_threshold']) / 2 - ) * bin_precsnow[:, t_start:t_stop] - self.bin_acc[:, t_start:t_stop] = bin_precsnow[:, t_start:t_stop] - self.bin_prec[:, t_start:t_stop] + self.bin_prec[:, t_start : t_stop + 1] = ( + 0.5 + (self.bin_temp[:, t_start : t_stop + 1] - self.modelprms['tsnow_threshold']) / 2 + ) * bin_precsnow[:, t_start : t_stop + 1] + self.bin_acc[:, t_start : t_stop + 1] = ( + bin_precsnow[:, t_start : t_stop + 1] - self.bin_prec[:, t_start : t_stop + 1] + ) # if temperature above maximum threshold, then all rain ( - self.bin_prec[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] > self.modelprms['tsnow_threshold'] + 1 + self.bin_prec[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] + 1 ] - ) = bin_precsnow[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] > self.modelprms['tsnow_threshold'] + 1 + ) = bin_precsnow[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] + 1 ] ( - self.bin_acc[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] > self.modelprms['tsnow_threshold'] + 1 + self.bin_acc[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] > self.modelprms['tsnow_threshold'] + 1 ] ) = 0 # if temperature below minimum threshold, then all snow ( - self.bin_acc[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] <= self.modelprms['tsnow_threshold'] - 1 + self.bin_acc[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] - 1 ] - ) = bin_precsnow[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] <= self.modelprms['tsnow_threshold'] - 1 + ) = bin_precsnow[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] - 1 ] ( - self.bin_prec[:, t_start:t_stop][ - self.bin_temp[:, t_start:t_stop] <= self.modelprms['tsnow_threshold'] - 1 + self.bin_prec[:, t_start : t_stop + 1][ + self.bin_temp[:, t_start : t_stop + 1] <= self.modelprms['tsnow_threshold'] - 1 ] ) = 0 @@ -668,8 +678,8 @@ def get_annual_mb( # calculate annually and place potential refreeze in user defined month if step == t_start: bin_temp_annual = annualweightedmean_array( - self.bin_temp[:, t_start:t_stop], - self.dates_table.iloc[t_start:t_stop, :], + self.bin_temp[:, t_start : t_stop + 1], + self.dates_table.iloc[t_start : t_stop + 1, :], ) bin_refreezepotential_annual = (-0.69 * bin_temp_annual + 0.0096) / 100 # Remove negative refreezing values @@ -751,7 +761,7 @@ def get_annual_mb( # ===== RETURN TO ANNUAL LOOP ===== # SURFACE TYPE (-) # Annual climatic mass balance [m w.e.] used to determine the surface type - self.glac_bin_massbalclim_annual[:, year_idx] = self.glac_bin_massbalclim[:, t_start:t_stop].sum(1) + self.glac_bin_massbalclim_annual[:, year_idx] = self.glac_bin_massbalclim[:, t_start : t_stop + 1].sum(1) # Update surface type for each bin self.surfacetype, firnline_idx = self._surfacetypebinsannual( self.surfacetype, self.glac_bin_massbalclim_annual, year_idx @@ -770,10 +780,10 @@ def get_annual_mb( option_areaconstant=option_areaconstant, ) - # Mass balance for each bin [m ice per second] - seconds_in_year = self.days_in_step[t_start:t_stop].sum() * 24 * 3600 + # Calculate mass balance rate for each bin (convert from [m w.e. yr^-1] to [m ice s^-1]) + seconds_in_year = self.days_in_step[t_start : t_stop + 1].sum() * 24 * 3600 mb = ( - self.glac_bin_massbalclim[:, t_start:t_stop].sum(1) + self.glac_bin_massbalclim[:, t_start : t_stop + 1].sum(1) * pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice'] / seconds_in_year @@ -829,7 +839,8 @@ def _convert_glacwide_results( """ # Glacier area glac_idx = glacier_area.nonzero()[0] - glacier_area_steps = glacier_area[:, np.newaxis].repeat(t_stop - t_start, axis=1) + # repeat glacier area for all time steps in year + glacier_area_steps = glacier_area[:, np.newaxis].repeat((t_stop + 1) - t_start, axis=1) # Check if need to adjust for complete removal of the glacier # - needed for accurate runoff calcs and accurate mass balance components @@ -846,7 +857,7 @@ def _convert_glacwide_results( ) # Check annual climatic mass balance (mwea) mb_mwea = ( - glacier_area * self.glac_bin_massbalclim[:, t_start:t_stop].sum(1) + glacier_area * self.glac_bin_massbalclim[:, t_start : t_stop + 1].sum(1) ).sum() / glacier_area.sum() else: mb_max_loss = 0 @@ -867,65 +878,65 @@ def _convert_glacwide_results( # Glacier-wide area (m2) self.glac_wide_area_annual[year_idx] = glacier_area.sum() # Glacier-wide temperature (degC) - self.glac_wide_temp[t_start:t_stop] = ( - self.bin_temp[:, t_start:t_stop][glac_idx] * glacier_area_steps[glac_idx] + self.glac_wide_temp[t_start : t_stop + 1] = ( + self.bin_temp[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) / glacier_area.sum() # Glacier-wide precipitation (m3) - self.glac_wide_prec[t_start:t_stop] = ( - self.bin_prec[:, t_start:t_stop][glac_idx] * glacier_area_steps[glac_idx] + self.glac_wide_prec[t_start : t_stop + 1] = ( + self.bin_prec[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) # Glacier-wide accumulation (m3 w.e.) - self.glac_wide_acc[t_start:t_stop] = ( - self.bin_acc[:, t_start:t_stop][glac_idx] * glacier_area_steps[glac_idx] + self.glac_wide_acc[t_start : t_stop + 1] = ( + self.bin_acc[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) # Glacier-wide refreeze (m3 w.e.) - self.glac_wide_refreeze[t_start:t_stop] = ( - self.glac_bin_refreeze[:, t_start:t_stop][glac_idx] * glacier_area_steps[glac_idx] + self.glac_wide_refreeze[t_start : t_stop + 1] = ( + self.glac_bin_refreeze[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) # Glacier-wide melt (m3 w.e.) - self.glac_wide_melt[t_start:t_stop] = ( - self.glac_bin_melt[:, t_start:t_stop][glac_idx] * glacier_area_steps[glac_idx] + self.glac_wide_melt[t_start : t_stop + 1] = ( + self.glac_bin_melt[:, t_start : t_stop + 1][glac_idx] * glacier_area_steps[glac_idx] ).sum(0) # Glacier-wide total mass balance (m3 w.e.) - self.glac_wide_massbaltotal[t_start:t_stop] = ( - self.glac_wide_acc[t_start:t_stop] - + self.glac_wide_refreeze[t_start:t_stop] - - self.glac_wide_melt[t_start:t_stop] - - self.glac_wide_frontalablation[t_start:t_stop] + self.glac_wide_massbaltotal[t_start : t_stop + 1] = ( + self.glac_wide_acc[t_start : t_stop + 1] + + self.glac_wide_refreeze[t_start : t_stop + 1] + - self.glac_wide_melt[t_start : t_stop + 1] + - self.glac_wide_frontalablation[t_start : t_stop + 1] ) # If mass loss more negative than glacier mass, reduce melt so glacier completely melts (no excess) if icethickness_t0 is not None and mb_mwea < mb_max_loss: - melt_yr_raw = self.glac_wide_melt[t_start:t_stop].sum() + melt_yr_raw = self.glac_wide_melt[t_start : t_stop + 1].sum() melt_yr_max = ( self.glac_wide_volume_annual[year_idx] * pygem_prms['constants']['density_ice'] / pygem_prms['constants']['density_water'] - + self.glac_wide_acc[t_start:t_stop].sum() - + self.glac_wide_refreeze[t_start:t_stop].sum() + + self.glac_wide_acc[t_start : t_stop + 1].sum() + + self.glac_wide_refreeze[t_start : t_stop + 1].sum() ) melt_frac = melt_yr_max / melt_yr_raw # Update glacier-wide melt (m3 w.e.) - self.glac_wide_melt[t_start:t_stop] = self.glac_wide_melt[t_start:t_stop] * melt_frac + self.glac_wide_melt[t_start : t_stop + 1] = self.glac_wide_melt[t_start : t_stop + 1] * melt_frac # Glacier-wide runoff (m3) - self.glac_wide_runoff[t_start:t_stop] = ( - self.glac_wide_prec[t_start:t_stop] - + self.glac_wide_melt[t_start:t_stop] - - self.glac_wide_refreeze[t_start:t_stop] + self.glac_wide_runoff[t_start : t_stop + 1] = ( + self.glac_wide_prec[t_start : t_stop + 1] + + self.glac_wide_melt[t_start : t_stop + 1] + - self.glac_wide_refreeze[t_start : t_stop + 1] ) # Snow line altitude (m a.s.l.) - heights_steps = heights[:, np.newaxis].repeat(t_stop - t_start, axis=1) + heights_steps = heights[:, np.newaxis].repeat((t_stop + 1) - t_start, axis=1) snow_mask = np.zeros(heights_steps.shape) - snow_mask[self.glac_bin_snowpack[:, t_start:t_stop] > 0] = 1 + snow_mask[self.glac_bin_snowpack[:, t_start : t_stop + 1] > 0] = 1 heights_steps_wsnow = heights_steps * snow_mask heights_steps_wsnow[heights_steps_wsnow == 0] = np.nan heights_change = np.zeros(heights.shape) heights_change[0:-1] = heights[0:-1] - heights[1:] try: snowline_idx = np.nanargmin(heights_steps_wsnow, axis=0) - self.glac_wide_snowline[t_start:t_stop] = heights[snowline_idx] - heights_change[snowline_idx] / 2 + self.glac_wide_snowline[t_start : t_stop + 1] = heights[snowline_idx] - heights_change[snowline_idx] / 2 except: snowline_idx = np.zeros((heights_steps_wsnow.shape[1])).astype(int) snowline_idx_nan = [] @@ -937,7 +948,7 @@ def _convert_glacwide_results( heights_manual = heights[snowline_idx] - heights_change[snowline_idx] / 2 heights_manual[snowline_idx_nan] = np.nan # this line below causes a potential All-NaN slice encountered issue at some time steps - self.glac_wide_snowline[t_start:t_stop] = heights_manual + self.glac_wide_snowline[t_start : t_stop + 1] = heights_manual # Equilibrium line altitude (m a.s.l.) ela_mask = np.zeros(heights.shape) @@ -954,33 +965,33 @@ def _convert_glacwide_results( offglac_idx = np.where(self.offglac_bin_area_annual[:, year_idx] > 0)[0] if option_areaconstant == False and len(offglac_idx) > 0: offglacier_area_steps = self.offglac_bin_area_annual[:, year_idx][:, np.newaxis].repeat( - t_stop - t_start, axis=1 + (t_stop + 1) - t_start, axis=1 ) # Off-glacier precipitation (m3) - self.offglac_wide_prec[t_start:t_stop] = ( - self.bin_prec[:, t_start:t_stop][offglac_idx] * offglacier_area_steps[offglac_idx] + self.offglac_wide_prec[t_start : t_stop + 1] = ( + self.bin_prec[:, t_start : t_stop + 1][offglac_idx] * offglacier_area_steps[offglac_idx] ).sum(0) # Off-glacier melt (m3 w.e.) - self.offglac_wide_melt[t_start:t_stop] = ( - self.offglac_bin_melt[:, t_start:t_stop][offglac_idx] * offglacier_area_steps[offglac_idx] + self.offglac_wide_melt[t_start : t_stop + 1] = ( + self.offglac_bin_melt[:, t_start : t_stop + 1][offglac_idx] * offglacier_area_steps[offglac_idx] ).sum(0) # Off-glacier refreeze (m3 w.e.) - self.offglac_wide_refreeze[t_start:t_stop] = ( - self.offglac_bin_refreeze[:, t_start:t_stop][offglac_idx] * offglacier_area_steps[offglac_idx] + self.offglac_wide_refreeze[t_start : t_stop + 1] = ( + self.offglac_bin_refreeze[:, t_start : t_stop + 1][offglac_idx] * offglacier_area_steps[offglac_idx] ).sum(0) # Off-glacier runoff (m3) - self.offglac_wide_runoff[t_start:t_stop] = ( - self.offglac_wide_prec[t_start:t_stop] - + self.offglac_wide_melt[t_start:t_stop] - - self.offglac_wide_refreeze[t_start:t_stop] + self.offglac_wide_runoff[t_start : t_stop + 1] = ( + self.offglac_wide_prec[t_start : t_stop + 1] + + self.offglac_wide_melt[t_start : t_stop + 1] + - self.offglac_wide_refreeze[t_start : t_stop + 1] ) # Off-glacier snowpack (m3 w.e.) - self.offglac_wide_snowpack[t_start:t_stop] = ( - self.offglac_bin_snowpack[:, t_start:t_stop][offglac_idx] * offglacier_area_steps[offglac_idx] + self.offglac_wide_snowpack[t_start : t_stop + 1] = ( + self.offglac_bin_snowpack[:, t_start : t_stop + 1][offglac_idx] * offglacier_area_steps[offglac_idx] ).sum(0) - def ensure_mass_conservation(self, diag, dates_table): + def ensure_mass_conservation(self, diag): """ 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 @@ -994,24 +1005,23 @@ def ensure_mass_conservation(self, diag, dates_table): total volume change and therefore do not impose limitations like this because they do not estimate the flux divergence. As a result, they may systematically overestimate mass loss compared to OGGM's dynamical model. """ - years = list(np.unique(dates_table['year'])) + years = list(np.unique(self.dates_table['year'])) # Compute annual volume change and melt values needed for adjustments vol_change_annual_mbmod = np.zeros(len(years)) vol_change_annual_mbmod_melt = np.zeros(len(years)) for nyear, year in enumerate(years): - step_idxs = np.where(self.dates_table.year == int(year)) - t_start = step_idxs[0][0] - t_stop = step_idxs[0][-1] + 1 + # get time step indices - note indexing should be [t_start:t_stop+1] to include final step in year + t_start, t_stop = self.get_step_inds(year) vol_change_annual_mbmod[nyear] = ( - self.glac_wide_massbaltotal[t_start:t_stop].sum() + self.glac_wide_massbaltotal[t_start : t_stop + 1].sum() * pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice'] ) vol_change_annual_mbmod_melt[nyear] = ( - self.glac_wide_melt[t_start:t_stop].sum() + self.glac_wide_melt[t_start : t_stop + 1].sum() * pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice'] ) @@ -1033,13 +1043,11 @@ def ensure_mass_conservation(self, diag, dates_table): 1 - vol_change_annual_dif[chg_idx_melt] / vol_change_annual_mbmod_melt[chg_idx_melt] ) - vol_change_annual_melt_reduction_monthly = np.zeros(dates_table.shape[0]) + vol_change_annual_melt_reduction_monthly = np.zeros(self.dates_table.shape[0]) for nyear, year in enumerate(years): - step_idxs = np.where(self.dates_table.year == int(year)) - t_start = step_idxs[0][0] - t_stop = step_idxs[0][-1] + 1 - - vol_change_annual_melt_reduction_monthly[t_start:t_stop] = vol_change_annual_melt_reduction[nyear] + # get time step indices - note indexing should be [t_start:t_stop+1] to include final step in year + t_start, t_stop = self.get_step_inds(year) + vol_change_annual_melt_reduction_monthly[t_start : t_stop + 1] = vol_change_annual_melt_reduction[nyear] # Glacier-wide melt (m3 w.e.) self.glac_wide_melt = self.glac_wide_melt * vol_change_annual_melt_reduction_monthly diff --git a/pygem/plot/graphics.py b/pygem/plot/graphics.py index 74f0bbba..103725b3 100644 --- a/pygem/plot/graphics.py +++ b/pygem/plot/graphics.py @@ -18,7 +18,15 @@ from pygem.utils.stats import effective_n -def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): +def plot_modeloutput_section( + model=None, + ax=None, + title='', + lnlabel=None, + legendon=True, + lgdkwargs={'loc': 'upper right', 'fancybox': False, 'borderaxespad': 0, 'handlelength': 1}, + **kwargs, +): """Plots the result of the model output along the flowline. A paired down version of OGGMs graphics.plot_modeloutput_section() @@ -40,15 +48,12 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): ax = fig.add_axes([0.07, 0.08, 0.7, 0.84]) else: fig = plt.gcf() + # get n lines plotted on figure + nlines = len(plt.gca().get_lines()) - # Compute area histo - area = np.array([]) height = np.array([]) bed = np.array([]) for cls in fls: - a = cls.widths_m * cls.dx_meter * 1e-6 - a = np.where(cls.thick > 0, a, 0) - area = np.concatenate((area, a)) height = np.concatenate((height, cls.surface_h)) bed = np.concatenate((bed, cls.bed_h)) ylim = [bed.min(), height.max()] @@ -57,8 +62,11 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): cls = fls[-1] x = np.arange(cls.nx) * cls.dx * cls.map_dx - # Plot the bed - ax.plot(x, cls.bed_h, color='k', linewidth=2.5, label='Bed (Parab.)') + if nlines == 0: + if getattr(model, 'do_calving', False): + ax.hlines(model.water_level, x[0], x[-1], linestyles=':', label='Water level', color='C0') + # Plot the bed + ax.plot(x, cls.bed_h, color='k', linewidth=2.5, label='Bed (Parab.)') # Plot glacier t1 = cls.thick[:-2] @@ -77,7 +85,7 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): else: srfls = '-' - ax.plot(x, cls.surface_h, color=srfcolor, linewidth=2, ls=srfls, label='Glacier') + ax.plot(x, cls.surface_h, color=srfcolor, linewidth=2, ls=srfls, label=lnlabel) # Plot tributaries for i, inflow in zip(cls.inflow_indices, cls.inflows): @@ -99,16 +107,14 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs): markeredgecolor='k', label='Tributary (inactive)', ) - if getattr(model, 'do_calving', False): - ax.hlines(model.water_level, x[0], x[-1], linestyles=':', color='C0') ax.set_ylim(ylim) - ax.spines['top'].set_color('none') ax.xaxis.set_ticks_position('bottom') ax.set_xlabel('Distance along flowline (m)') ax.set_ylabel('Altitude (m)') - + if legendon: + ax.legend(**lgdkwargs) # Title ax.set_title(title, loc='left') diff --git a/pygem/pygem_modelsetup.py b/pygem/pygem_modelsetup.py index ea215f9f..0edeb9a4 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -135,52 +135,6 @@ def datesmodelrun( return dates_table -def daysinmonth(year, month): - """ - Return days in month based on the month and year - - Parameters - ---------- - year : str - month : str - - Returns - ------- - integer of the days in the month - """ - if year % 4 == 0: - daysinmonth_dict = { - 1: 31, - 2: 29, - 3: 31, - 4: 30, - 5: 31, - 6: 30, - 7: 31, - 8: 31, - 9: 30, - 10: 31, - 11: 30, - 12: 31, - } - else: - daysinmonth_dict = { - 1: 31, - 2: 28, - 3: 31, - 4: 30, - 5: 31, - 6: 30, - 7: 31, - 8: 31, - 9: 30, - 10: 31, - 11: 30, - 12: 31, - } - return daysinmonth_dict[month] - - def hypsometrystats(hyps_table, thickness_table): """Calculate the volume and mean associated with the hypsometry data. diff --git a/pygem/tests/test_04_auxiliary.py b/pygem/tests/test_04_auxiliary.py new file mode 100644 index 00000000..c9db58fc --- /dev/null +++ b/pygem/tests/test_04_auxiliary.py @@ -0,0 +1,87 @@ +import glob +import os +import subprocess + +import pytest + +from pygem.setup.config import ConfigManager + +""" +Test suite to any necessary aux""" + + +@pytest.fixture(scope='module') +def rootdir(): + config_manager = ConfigManager() + pygem_prms = config_manager.read_config() + return pygem_prms['root'] + + +def test_simulation_massredistribution_dynamics(rootdir): + """ + Test the run_simulation CLI script with the "MassRedistributionCurves" dynamical option. + """ + + # Run run_simulation CLI script + subprocess.run( + [ + 'run_simulation', + '-rgi_glac_number', + '1.03622', + '-option_calibration', + 'MCMC', + '-sim_climate_name', + 'ERA5', + '-sim_startyear', + '2000', + '-sim_endyear', + '2019', + '-nsims', + '1', + '-option_dynamics', + 'MassRedistributionCurves', + '-outputfn_sfix', + 'mrcdynamics_', + ], + check=True, + ) + + # Check if output files were created + outdir = os.path.join(rootdir, 'Output', 'simulations', '01', 'ERA5') + output_files = glob.glob(os.path.join(outdir, '**', '*_mrcdynamics_all.nc'), recursive=True) + assert output_files, f'Simulation output file not found in {outdir}' + + +def test_simulation_no_dynamics(rootdir): + """ + Test the run_simulation CLI script with no dynamics option. + """ + + # Run run_simulation CLI script + subprocess.run( + [ + 'run_simulation', + '-rgi_glac_number', + '1.03622', + '-option_calibration', + 'MCMC', + '-sim_climate_name', + 'ERA5', + '-sim_startyear', + '2000', + '-sim_endyear', + '2019', + '-nsims', + '1', + '-option_dynamics', + 'None', + '-outputfn_sfix', + 'nodynamics_', + ], + check=True, + ) + + # Check if output files were created + outdir = os.path.join(rootdir, 'Output', 'simulations', '01', 'ERA5') + output_files = glob.glob(os.path.join(outdir, '**', '*_nodynamics_all.nc'), recursive=True) + assert output_files, f'Simulation output file not found in {outdir}' diff --git a/pygem/tests/test_04_postproc.py b/pygem/tests/test_05_postproc.py similarity index 96% rename from pygem/tests/test_04_postproc.py rename to pygem/tests/test_05_postproc.py index 297146a3..35a60eb8 100644 --- a/pygem/tests/test_04_postproc.py +++ b/pygem/tests/test_05_postproc.py @@ -99,7 +99,6 @@ def test_check_compiled_product(rootdir): vars_to_check = [item for item in vars_to_check if item not in vars_to_skip] for var in vars_to_check: - print(var) # skip mad if 'mad' in var: continue @@ -115,13 +114,13 @@ def test_check_compiled_product(rootdir): # pull data values simvals = simvar.values - compvals = compvar.values[0, :, :] # first index is the glacier index + compvals = compvar.values[0, :, :] # 0th index is the glacier index # check that compiled product has same shape as original data assert simvals.shape == compvals.shape, ( f'Compiled product shape {compvals.shape} does not match original data shape {simvals.shape}' ) # check that compiled product matches original data - assert np.allclose(simvals, compvals, rtol=1e-8, atol=1e-12), ( + assert np.allclose(simvals, compvals, rtol=1e-8, atol=1e-12, equal_nan=True), ( f'Compiled product for {var} does not match original data' )