From 69064137a366dd3e73adac96eb8566c3ea82cd45 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 15 Aug 2023 17:31:56 +0200 Subject: [PATCH 1/9] add option for monthly interpolation to distribute_thickness_from_simulation --- oggm/sandbox/distribute_2d.py | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index 6398db749..ba5a03f47 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -195,7 +195,8 @@ def assign_points_to_band(gdir, topo_variable='glacier_topo_smoothed', @entity_task(log, writes=['gridded_data']) def distribute_thickness_from_simulation(gdir, input_filesuffix='', ys=None, ye=None, - smooth_radius=None): + smooth_radius=None, + add_monthly=False): """Redistributes the simulated flowline area and volume back onto the 2D grid. For this to work, the glacier cannot advance beyond its initial area! @@ -230,6 +231,9 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', pixel size of the gaussian smoothing. Default is to use cfg.PARAMS['smooth_window'] (i.e. a size in meters). Set to zero to suppress smoothing. + add_monthly : bool + If True the yearly flowline diagnostics will be linearly interpolated + to a monthly resolution. Default is False """ fp = gdir.get_filepath('fl_diagnostics', filesuffix=input_filesuffix) @@ -240,6 +244,20 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', dg = dg.sel(time=slice(ye, ye)) dg = dg.load() + if add_monthly: + # create new monthly time coordinate, last year only with month 1 + years = np.append(np.repeat(dg.time[:-1], 12), + dg.time[-1]) + months = np.append(np.tile(np.arange(1, 13), len(dg.time[:-1])), + 1) + time_monthly = utils.date_to_floatyear(years, months) + + # interpolate and add years and months as new coords + dg = dg[['area_m2', 'volume_m3']].interp(time=time_monthly, + method='linear') + dg.coords['calender_year'] = ('time', years) + dg.coords['calender_month'] = ('time', months) + with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: band_index_mask = ds.band_index.data rank_per_band = ds.rank_per_band.data @@ -292,9 +310,19 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: ds = ds.load() - ds['time'] = dg['time'] + # the distinction between time_monthly and time is needed to have data in + # yearly AND monthly resolution in the same gridded data, maybe and one + # point we decide for one option + if add_monthly: + time_var = 'time_monthly' + else: + time_var = 'time' + ds.coords[time_var] = dg['time'].data vn = "simulation_distributed_thickness" + input_filesuffix if vn in ds: warnings.warn(f'Overwriting existing variable {vn}') - ds[vn] = (('time', 'y', 'x',), out_thick) + ds[vn] = ((time_var, 'y', 'x',), out_thick) + if add_monthly: + ds.coords['calender_year_monthly'] = (time_var, dg.calender_year.data) + ds.coords['calender_month_monthly'] = (time_var, dg.calender_month.data) ds.to_netcdf(gdir.get_filepath('gridded_data')) From cc65d674a58ec59810ab49a979d9ba601f5ed132 Mon Sep 17 00:00:00 2001 From: afisc Date: Wed, 16 Aug 2023 14:57:35 +0200 Subject: [PATCH 2/9] added flowline thickness thresholding to the distribute_2d functionality --- oggm/sandbox/distribute_2d.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index ba5a03f47..417f0ba8f 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -196,7 +196,8 @@ def assign_points_to_band(gdir, topo_variable='glacier_topo_smoothed', def distribute_thickness_from_simulation(gdir, input_filesuffix='', ys=None, ye=None, smooth_radius=None, - add_monthly=False): + add_monthly=False, + fl_thickness_threshold=0): """Redistributes the simulated flowline area and volume back onto the 2D grid. For this to work, the glacier cannot advance beyond its initial area! @@ -234,6 +235,9 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', add_monthly : bool If True the yearly flowline diagnostics will be linearly interpolated to a monthly resolution. Default is False + fl_thickness_threshold: int + A minimum threshold (all values below the threshold are set to 0) is applied to the area and volume of the + flowline diagnostics before the distribution process. """ fp = gdir.get_filepath('fl_diagnostics', filesuffix=input_filesuffix) @@ -244,6 +248,11 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', dg = dg.sel(time=slice(ye, ye)) dg = dg.load() + # tackling flickering issue + if fl_thickness_threshold: + dg['area_m2'].values = np.where(dg['thickness_m'] < fl_thickness_threshold, 0, dg['area_m2']) + dg['volume_m3'].values = np.where(dg['thickness_m'] < fl_thickness_threshold, 0, dg['volume_m3']) + if add_monthly: # create new monthly time coordinate, last year only with month 1 years = np.append(np.repeat(dg.time[:-1], 12), From 9867741c8976e6c83943e716684c2e00c838ffaf Mon Sep 17 00:00:00 2001 From: afisc Date: Thu, 17 Aug 2023 14:00:54 +0200 Subject: [PATCH 3/9] added area smoothing for a smoother animation --- oggm/sandbox/distribute_2d.py | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index 417f0ba8f..076a5b93c 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -197,7 +197,8 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', ys=None, ye=None, smooth_radius=None, add_monthly=False, - fl_thickness_threshold=0): + fl_thickness_threshold=0, + area_smoothing = False): """Redistributes the simulated flowline area and volume back onto the 2D grid. For this to work, the glacier cannot advance beyond its initial area! @@ -238,6 +239,11 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', fl_thickness_threshold: int A minimum threshold (all values below the threshold are set to 0) is applied to the area and volume of the flowline diagnostics before the distribution process. + area_smoothing: bool + If True, the area and volume of the flowline diagnostics will be smoothed linearly for each grid point of the + flowline. It is simply a linear interpolation between the beginning area/volume of the grid point and its + state of complete mass loss. This is not area and volume conserving anymore and can also not represent regain + of glacier mass during the simulation. """ fp = gdir.get_filepath('fl_diagnostics', filesuffix=input_filesuffix) @@ -248,11 +254,26 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', dg = dg.sel(time=slice(ye, ye)) dg = dg.load() + # smoothing area + + # tackling flickering issue if fl_thickness_threshold: dg['area_m2'].values = np.where(dg['thickness_m'] < fl_thickness_threshold, 0, dg['area_m2']) dg['volume_m3'].values = np.where(dg['thickness_m'] < fl_thickness_threshold, 0, dg['volume_m3']) + + if area_smoothing: + for variable_name in ['area_m2', 'volume_m3']: + for i in np.arange(0, len(dg.dis_along_flowline)): + nr_non_zero = sum(dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i]}] != 0).item() + 1 + if nr_non_zero > len(dg.time): + nr_non_zero -= 1 + new_values = np.linspace(dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i]}][0], 0, + num=nr_non_zero) + dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i], + 'time': slice(dg.time[nr_non_zero - 1])}] = new_values + if add_monthly: # create new monthly time coordinate, last year only with month 1 years = np.append(np.repeat(dg.time[:-1], 12), From b7fac85b2e887fcdf1e707ca7fb6b509354591c1 Mon Sep 17 00:00:00 2001 From: afisc Date: Thu, 17 Aug 2023 17:42:45 +0200 Subject: [PATCH 4/9] added rolling_window option to area_smoothing --- oggm/sandbox/distribute_2d.py | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index 076a5b93c..d2f0c3d02 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -198,7 +198,8 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', smooth_radius=None, add_monthly=False, fl_thickness_threshold=0, - area_smoothing = False): + area_smoothing = False, + rolling_window = 0): """Redistributes the simulated flowline area and volume back onto the 2D grid. For this to work, the glacier cannot advance beyond its initial area! @@ -244,6 +245,10 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', flowline. It is simply a linear interpolation between the beginning area/volume of the grid point and its state of complete mass loss. This is not area and volume conserving anymore and can also not represent regain of glacier mass during the simulation. + rolling_window: int + 'area_smoothing' has to be switched on for this arg to have an effect. + If a rolling window value is given, the area smoothing is done with the + rolling window algorithm instead of linear interpolation. """ fp = gdir.get_filepath('fl_diagnostics', filesuffix=input_filesuffix) @@ -254,8 +259,6 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', dg = dg.sel(time=slice(ye, ye)) dg = dg.load() - # smoothing area - # tackling flickering issue if fl_thickness_threshold: @@ -264,15 +267,19 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', if area_smoothing: - for variable_name in ['area_m2', 'volume_m3']: - for i in np.arange(0, len(dg.dis_along_flowline)): - nr_non_zero = sum(dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i]}] != 0).item() + 1 - if nr_non_zero > len(dg.time): - nr_non_zero -= 1 - new_values = np.linspace(dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i]}][0], 0, - num=nr_non_zero) - dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i], - 'time': slice(dg.time[nr_non_zero - 1])}] = new_values + if rolling_window: + dg[['area_m2', 'volume_m3']] = dg[['area_m2', 'volume_m3']].rolling(min_periods=1, time=rolling_window, + center=True).mean() + else: + for variable_name in ['area_m2', 'volume_m3']: + for i in np.arange(0, len(dg.dis_along_flowline)): + nr_non_zero = sum(dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i]}] != 0).item() + 1 + if nr_non_zero > len(dg.time): + nr_non_zero -= 1 + new_values = np.linspace(dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i]}][0], 0, + num=nr_non_zero) + dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i], + 'time': slice(dg.time[nr_non_zero - 1])}] = new_values if add_monthly: # create new monthly time coordinate, last year only with month 1 From b887d030347c3bdde8258572507a1bd501893f66 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 18 Aug 2023 11:12:28 +0200 Subject: [PATCH 5/9] code cleanup, added option for area evolution plot --- oggm/sandbox/distribute_2d.py | 76 +++++++++++++++++++---------------- 1 file changed, 41 insertions(+), 35 deletions(-) diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index d2f0c3d02..1e6c0dc63 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -9,6 +9,7 @@ from scipy.stats import mstats from oggm.core.gis import gaussian_blur from oggm.utils import ncDataset, entity_task +import matplotlib.pyplot as plt # Module logger log = logging.getLogger(__name__) @@ -198,8 +199,8 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', smooth_radius=None, add_monthly=False, fl_thickness_threshold=0, - area_smoothing = False, - rolling_window = 0): + rolling_mean_smoothing=False, + show_area_plot=False): """Redistributes the simulated flowline area and volume back onto the 2D grid. For this to work, the glacier cannot advance beyond its initial area! @@ -237,18 +238,20 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', add_monthly : bool If True the yearly flowline diagnostics will be linearly interpolated to a monthly resolution. Default is False - fl_thickness_threshold: int - A minimum threshold (all values below the threshold are set to 0) is applied to the area and volume of the - flowline diagnostics before the distribution process. - area_smoothing: bool - If True, the area and volume of the flowline diagnostics will be smoothed linearly for each grid point of the - flowline. It is simply a linear interpolation between the beginning area/volume of the grid point and its - state of complete mass loss. This is not area and volume conserving anymore and can also not represent regain - of glacier mass during the simulation. - rolling_window: int - 'area_smoothing' has to be switched on for this arg to have an effect. - If a rolling window value is given, the area smoothing is done with the - rolling window algorithm instead of linear interpolation. + fl_thickness_threshold : float + A minimum threshold (all values below the threshold are set to 0) is + applied to the area and volume of the flowline diagnostics before the + distribution process. Default is 0 (using no threshold). + rolling_mean_smoothing : bool or int + If int, the area and volume of the flowline diagnostics will be + smoothed using a rolling mean. The window size is defined with this + number. If True a default window size of 3 is used. If False no + smoothing is applied. Default is False. + show_area_plot : bool + If True, only a plot of the 'original-total-area- evolution' and the + 'smoothed-total-area- evolution' (after using fl_thickness_threshold + and rolling_mean_smoothing) is returned. This is useful for finding the + best smoothing parameters for the visualisation of your glacier. """ fp = gdir.get_filepath('fl_diagnostics', filesuffix=input_filesuffix) @@ -259,28 +262,22 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', dg = dg.sel(time=slice(ye, ye)) dg = dg.load() + # save the original area evolution for the area plot + if show_area_plot: + area_evolution_orig = dg['area_m2'].sum(dim='dis_along_flowline') - # tackling flickering issue - if fl_thickness_threshold: - dg['area_m2'].values = np.where(dg['thickness_m'] < fl_thickness_threshold, 0, dg['area_m2']) - dg['volume_m3'].values = np.where(dg['thickness_m'] < fl_thickness_threshold, 0, dg['volume_m3']) + # applying the thickness threshold + dg = xr.where(dg['thickness_m'] < fl_thickness_threshold, 0, dg) + # applying the rolling mean smoothing + if rolling_mean_smoothing: + if isinstance(rolling_mean_smoothing, bool): + rolling_mean_smoothing = 3 - if area_smoothing: - if rolling_window: - dg[['area_m2', 'volume_m3']] = dg[['area_m2', 'volume_m3']].rolling(min_periods=1, time=rolling_window, - center=True).mean() - else: - for variable_name in ['area_m2', 'volume_m3']: - for i in np.arange(0, len(dg.dis_along_flowline)): - nr_non_zero = sum(dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i]}] != 0).item() + 1 - if nr_non_zero > len(dg.time): - nr_non_zero -= 1 - new_values = np.linspace(dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i]}][0], 0, - num=nr_non_zero) - dg[variable_name].loc[{'dis_along_flowline': dg.dis_along_flowline[i], - 'time': slice(dg.time[nr_non_zero - 1])}] = new_values + dg[['area_m2', 'volume_m3']] = dg[['area_m2', 'volume_m3']].rolling( + min_periods=1, time=rolling_mean_smoothing, center=True).mean() + # monthly interpolation for higher temporal resolution if add_monthly: # create new monthly time coordinate, last year only with month 1 years = np.append(np.repeat(dg.time[:-1], 12), @@ -295,13 +292,21 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', dg.coords['calender_year'] = ('time', years) dg.coords['calender_month'] = ('time', months) + if show_area_plot: + area_evolution_orig.plot(label='original') + dg['area_m2'].sum(dim='dis_along_flowline').plot(label='smoothed') + plt.legend() + plt.show() + return None + with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: band_index_mask = ds.band_index.data rank_per_band = ds.rank_per_band.data glacier_mask = ds.glacier_mask.data == 1 orig_distrib_thick = ds.distributed_thickness.data - band_ids, counts = np.unique(np.sort(band_index_mask[glacier_mask]), return_counts=True) + band_ids, counts = np.unique(np.sort(band_index_mask[glacier_mask]), + return_counts=True) dx2 = gdir.grid.dx**2 out_thick = np.zeros((len(dg.time), *glacier_mask.shape)) @@ -318,7 +323,8 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', if band_area != 0: # We have some ice left pix_cov = (band_area / dx2) + residual_pix - mask = (band_index_mask == band_id) & (rank_per_band >= (npix - pix_cov)) + mask = (band_index_mask == band_id) & \ + (rank_per_band >= (npix - pix_cov)) residual_pix = pix_cov - mask.sum() vol_orig = np.where(mask, orig_distrib_thick, 0).sum() * dx2 area_dis = mask.sum() * dx2 @@ -348,7 +354,7 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', ds = ds.load() # the distinction between time_monthly and time is needed to have data in - # yearly AND monthly resolution in the same gridded data, maybe and one + # yearly AND monthly resolution in the same gridded data, maybe at one # point we decide for one option if add_monthly: time_var = 'time_monthly' From ed826a5ef63814ec3ddb00a6a58d73300f0905a7 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 18 Aug 2023 12:24:23 +0200 Subject: [PATCH 6/9] small change in redistribuion algorithm --- oggm/sandbox/distribute_2d.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index 1e6c0dc63..f1d7054ed 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -314,7 +314,6 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', dgy = dg.sel(time=yr) - residual_pix = 0 area_cov = 0 new_thick = out_thick[i, :] for band_id, npix in zip(band_ids.astype(int), counts): @@ -322,10 +321,9 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', band_volume = dgy.volume_m3.values[band_id] if band_area != 0: # We have some ice left - pix_cov = (band_area / dx2) + residual_pix + pix_cov = band_area / dx2 mask = (band_index_mask == band_id) & \ (rank_per_band >= (npix - pix_cov)) - residual_pix = pix_cov - mask.sum() vol_orig = np.where(mask, orig_distrib_thick, 0).sum() * dx2 area_dis = mask.sum() * dx2 thick_cor = (vol_orig - band_volume) / area_dis From 10822f129bf0bdb87276e9724ca0acd04c04b0ef Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 18 Aug 2023 17:00:32 +0200 Subject: [PATCH 7/9] added possibility to provide flowline diagnostics directly --- oggm/sandbox/distribute_2d.py | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index f1d7054ed..76bfd9e06 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -1,6 +1,7 @@ import logging import warnings +import oggm import oggm.cfg as cfg from oggm import utils import numpy as np @@ -195,6 +196,7 @@ def assign_points_to_band(gdir, topo_variable='glacier_topo_smoothed', @entity_task(log, writes=['gridded_data']) def distribute_thickness_from_simulation(gdir, input_filesuffix='', + fl_diag=None, ys=None, ye=None, smooth_radius=None, add_monthly=False, @@ -225,6 +227,10 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', where to write the data input_filesuffix : str the filesuffix of the flowline diagnostics file. + fl_diag : xarray.core.dataset.Dataset + can directly provide a flowline diagnostics file. If provided + 'input_filesuffix' is only used for the name to save the distributed + data in gridded data. ys : int pick another year to start the series (default: the first year of the diagnostic file) @@ -254,13 +260,16 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', best smoothing parameters for the visualisation of your glacier. """ - fp = gdir.get_filepath('fl_diagnostics', filesuffix=input_filesuffix) - with xr.open_dataset(fp) as dg: - assert len(dg.flowlines.data) == 1, 'Only works with one flowline.' - with xr.open_dataset(fp, group=f'fl_0') as dg: - if ys or ye: - dg = dg.sel(time=slice(ye, ye)) - dg = dg.load() + if fl_diag is not None: + dg = fl_diag + else: + fp = gdir.get_filepath('fl_diagnostics', filesuffix=input_filesuffix) + with xr.open_dataset(fp) as dg: + assert len(dg.flowlines.data) == 1, 'Only works with one flowline.' + with xr.open_dataset(fp, group=f'fl_0') as dg: + if ys or ye: + dg = dg.sel(time=slice(ye, ye)) + dg = dg.load() # save the original area evolution for the area plot if show_area_plot: From caaa2d3d29fedde83adb1ce08141a1fa0ec0cf4f Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 21 Aug 2023 11:33:58 +0200 Subject: [PATCH 8/9] try to make adding_points_to_band more robuste --- oggm/sandbox/distribute_2d.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index 76bfd9e06..5955ece04 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -143,6 +143,17 @@ def assign_points_to_band(gdir, topo_variable='glacier_topo_smoothed', npix_per_band = fl.bin_area_m2 / (gdir.grid.dx ** 2) nnpix_per_band_cumsum = np.around(npix_per_band[::-1].cumsum()[::-1]) + # check if their is a differenze between the flowline area and gridded area + pix_diff = int(nnpix_per_band_cumsum[0] - len(topo_data_flat)) + if pix_diff > 0: + # we have more fl pixels, for this we can adapt a little + # search for the largest differences between elevation bands and reduce + # area for one pixel + sorted_indices = np.argsort(np.abs(np.diff(nnpix_per_band_cumsum))) + npix_per_band[sorted_indices[-pix_diff:]] = \ + npix_per_band[sorted_indices[-pix_diff:]] - 1 + nnpix_per_band_cumsum = np.around(npix_per_band[::-1].cumsum()[::-1]) + rank_elev = mstats.rankdata(topo_data_flat) bins = nnpix_per_band_cumsum[nnpix_per_band_cumsum > 0].copy() bins[0] = len(topo_data_flat) + 1 From cece69c90594cef9b2f91027cc1b42a9e88d809b Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 21 Aug 2023 12:08:32 +0200 Subject: [PATCH 9/9] use np.close to check for zero area in redistribute --- oggm/sandbox/distribute_2d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index 5955ece04..b79b4599d 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -339,7 +339,7 @@ def distribute_thickness_from_simulation(gdir, input_filesuffix='', for band_id, npix in zip(band_ids.astype(int), counts): band_area = dgy.area_m2.values[band_id] band_volume = dgy.volume_m3.values[band_id] - if band_area != 0: + if ~np.isclose(band_area, 0): # We have some ice left pix_cov = band_area / dx2 mask = (band_index_mask == band_id) & \