From 3dd30ba3735782f6b739bfc76197f9c2203cb810 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 26 Sep 2019 11:37:42 +0100 Subject: [PATCH 1/2] Use faster clips everywhere --- oggm/core/centerlines.py | 18 +++++++------- oggm/core/climate.py | 8 +++---- oggm/core/flowline.py | 24 +++++++++---------- oggm/core/gcm_climate.py | 2 +- oggm/core/gis.py | 15 ++++++------ oggm/core/inversion.py | 12 +++++----- oggm/core/massbalance.py | 11 ++++----- oggm/core/sia2d.py | 8 ++++--- oggm/tests/test_utils.py | 51 ++++++++++++++++++++++++++++++++++++++++ oggm/utils/_funcs.py | 15 ++++++++++++ 10 files changed, 116 insertions(+), 48 deletions(-) diff --git a/oggm/core/centerlines.py b/oggm/core/centerlines.py index f39d5c99c..9f7d6a2f6 100644 --- a/oggm/core/centerlines.py +++ b/oggm/core/centerlines.py @@ -121,11 +121,11 @@ def set_flows_to(self, other, check_tail=True, last_point=False): ind_closest = np.argmin(np.abs(other.dis_on_line - prdis)).item() n = len(other.dis_on_line) if n >= 9: - ind_closest = np.clip(ind_closest, 4, n-5) + ind_closest = utils.clip_scalar(ind_closest, 4, n-5) elif n >= 7: - ind_closest = np.clip(ind_closest, 3, n-4) + ind_closest = utils.clip_scalar(ind_closest, 3, n-4) elif n >= 5: - ind_closest = np.clip(ind_closest, 2, n-3) + ind_closest = utils.clip_scalar(ind_closest, 2, n-3) p = shpg.Point(other.line.coords[int(ind_closest)]) self.flows_to_point = p elif last_point: @@ -708,7 +708,7 @@ def _get_centerlines_heads(gdir, ext_yx, zoutline, single_fl, # Size of the half window to use to look for local maximas maxorder = np.rint(cfg.PARAMS['localmax_window'] / gdir.grid.dx) - maxorder = np.clip(maxorder, 5., np.rint((len(zoutline) / 5.))) + maxorder = utils.clip_scalar(maxorder, 5., np.rint((len(zoutline) / 5.))) heads_idx = scipy.signal.argrelmax(zoutline, mode='wrap', order=maxorder.astype(np.int64)) if single_fl or len(heads_idx[0]) <= 1: @@ -731,7 +731,7 @@ def _get_centerlines_heads(gdir, ext_yx, zoutline, single_fl, # get radius of the buffer according to Kienholz eq. (1) radius = cfg.PARAMS['q1'] * geom['polygon_area'] + cfg.PARAMS['q2'] - radius = np.clip(radius, 0, cfg.PARAMS['rmax']) + radius = utils.clip_scalar(radius, 0, cfg.PARAMS['rmax']) radius /= gdir.grid.dx # in raster coordinates # Plus our criteria, quite useful to remove short lines: radius += cfg.PARAMS['flowline_junction_pix'] * cfg.PARAMS['flowline_dx'] @@ -1113,7 +1113,7 @@ def _parabolic_bed_from_topo(gdir, idl, interpolator): # We forbid super small shapes (important! This can lead to huge volumes) # Sometimes the parabola fits in flat areas are very good, implying very # flat parabolas. - bed_int = bed_int.clip(cfg.PARAMS['downstream_min_shape'], None) + bed_int = utils.clip_min(bed_int, cfg.PARAMS['downstream_min_shape']) # Smoothing bed_ma = pd.Series(bed_int) @@ -1804,7 +1804,7 @@ def catchment_width_correction(gdir): # Interpolate widths widths = utils.interp_nans(fl.widths) - widths = np.clip(widths, 0.1, None) + widths = utils.clip_min(widths, 0.1) # Get topo per catchment and per flowline point fhgt = fl.surface_h @@ -1815,11 +1815,11 @@ def catchment_width_correction(gdir): # Sometimes, the centerline does not reach as high as each pix on the # glacier. (e.g. RGI40-11.00006) - catch_h = np.clip(catch_h, None, maxh) + catch_h = utils.clip_max(catch_h, maxh) # Same for min if fl.flows_to is None: # We clip only for main flowline (this has reasons) - catch_h = np.clip(catch_h, minh, None) + catch_h = utils.clip_min(catch_h, minh) # Now decide on a binsize which ensures at least N element per bin bsize = cfg.PARAMS['base_binsize'] diff --git a/oggm/core/climate.py b/oggm/core/climate.py index c95415a67..d6b25fcf2 100644 --- a/oggm/core/climate.py +++ b/oggm/core/climate.py @@ -311,7 +311,7 @@ def process_cru_data(gdir): ts_pre.values, ts_pre_ano.values) # The last step might create negative values (unlikely). Clip them - ts_pre.values = ts_pre.values.clip(0) + ts_pre.values = utils.clip_min(ts_pre.values, 0) # done loc_hgt = loc_hgt[1, 1] @@ -449,7 +449,7 @@ def process_dummy_cru_file(gdir, sigma_temp=2, sigma_prcp=0.5, seed=None): loc_pre = xr.DataArray(loc_pre[:, 1, 1], dims=['month'], coords={'month': np.arange(1, 13)}) - ts_pre = (rng.randn(len(time)) * sigma_prcp + 1).clip(0) + ts_pre = utils.clip_min(rng.randn(len(time)) * sigma_prcp + 1, 0) ts_pre = xr.DataArray(ts_pre, dims=['time'], coords={'time': time}) @@ -680,7 +680,7 @@ def mb_climate_on_height(gdir, heights, *, time_range=None, year_range=None): grad_temp *= (heights.repeat(len(time)).reshape(grad_temp.shape) - ref_hgt) temp2d = np.atleast_2d(itemp).repeat(npix, 0) + grad_temp temp2dformelt = temp2d - temp_melt - temp2dformelt = np.clip(temp2dformelt, 0, None) + temp2dformelt = utils.clip_min(temp2dformelt, 0) # Compute solid precipitation from total precipitation prcpsol = np.atleast_2d(iprcp).repeat(npix, 0) fac = 1 - (temp2d - temp_all_solid) / (temp_all_liq - temp_all_solid) @@ -1104,7 +1104,7 @@ def local_t_star(gdir, *, ref_df=None, tstar=None, bias=None): # Clip it? if cfg.PARAMS['clip_mu_star']: - mustar = np.clip(mustar, 0, None) + mustar = utils.clip_min(mustar, 0) # If mu out of bounds, raise if not (cfg.PARAMS['min_mu_star'] <= mustar <= cfg.PARAMS['max_mu_star']): diff --git a/oggm/core/flowline.py b/oggm/core/flowline.py index 18a626610..f2c6fca07 100644 --- a/oggm/core/flowline.py +++ b/oggm/core/flowline.py @@ -62,7 +62,7 @@ def __init__(self, line=None, dx=1, map_dx=None, super(Flowline, self).__init__(line, dx, surface_h) - self._thick = (surface_h - bed_h).clip(0.) + self._thick = utils.clip_min(surface_h - bed_h, 0.) self.map_dx = map_dx self.dx_meter = map_dx * self.dx self.bed_h = bed_h @@ -80,7 +80,7 @@ def thick(self): @thick.setter def thick(self, value): - self._thick = value.clip(0) + self._thick = utils.clip_min(value, 0) @Centerline.surface_h.getter def surface_h(self): @@ -1094,7 +1094,7 @@ def step(self, dt): # Time step self.dt_warning = dt < min_dt - dt = np.clip(dt, min_dt, self.max_dt) + dt = utils.clip_scalar(dt, min_dt, self.max_dt) # A second loop for the mass exchange for fl_id, fl in enumerate(self.fls): @@ -1121,19 +1121,19 @@ def step(self, dt): trib_flux*dt/dx + mb) # Keep positive values only and store - fl.section = new_section.clip(0) + fl.section = utils.clip_min(new_section, 0) # Add the last flux to the tributary # this works because the lines are sorted in order if is_trib: # tr tuple: line_index, start, stop, gaussian_kernel - self.trib_flux[tr[0]][tr[1]:tr[2]] += (flx_stag[-1].clip(0) * - tr[3]) + self.trib_flux[tr[0]][tr[1]:tr[2]] += \ + utils.clip_min(flx_stag[-1], 0) * tr[3] elif self.is_tidewater: # -2 because the last flux is zero per construction # not sure at all if this is the way to go but mass # conservation is OK - self.calving_m3_since_y0 += flx_stag[-2].clip(0)*dt + self.calving_m3_since_y0 += utils.clip_min(flx_stag[-2], 0)*dt # Next step self.t += dt @@ -1222,7 +1222,7 @@ def step(self, dt): # there can't be more negative mb than there is section # this isn't an exact solution unfortunately # TODO: exact solution for mass conservation - mb = mb.clip(-sec) + mb = utils.clip_min(mb, -sec) self.total_mass += np.sum(mb * dx) @@ -1265,7 +1265,7 @@ def step(self, dt): # This is to guarantee a precise arrival on a specific date if asked min_dt = dt if dt < self.min_dt else self.min_dt - dt = np.clip(dt, min_dt, self.max_dt) + dt = utils.clip_scalar(dt, min_dt, self.max_dt) fl = self.fls[0] dx = fl.dx_meter @@ -1309,7 +1309,7 @@ def step(self, dt): NewIceThickness[-1] = thick[fl.nx-2] - fl.thick = NewIceThickness.clip(0) + fl.thick = utils.clip_min(NewIceThickness, 0) # Next step self.t += dt @@ -1368,7 +1368,7 @@ def step(self, dt): # Guarantee a precise arrival on a specific date if asked min_dt = dt if dt < self.min_dt else self.min_dt - dt = np.clip(dt, min_dt, self.max_dt) + dt = utils.clip_scalar(dt, min_dt, self.max_dt) with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=RuntimeWarning) @@ -1705,7 +1705,7 @@ def init_present_time_glacier(gdir): # Interpolate bed_shape = utils.interp_nans(np.append(bed_shape, dic_ds['bedshapes'][0])) - bed_shape = bed_shape[:-1].clip(min_shape) + bed_shape = utils.clip_min(bed_shape[:-1], min_shape) # Correct the section volume h = inv['thick'] diff --git a/oggm/core/gcm_climate.py b/oggm/core/gcm_climate.py index 938fd1320..903c697f6 100644 --- a/oggm/core/gcm_climate.py +++ b/oggm/core/gcm_climate.py @@ -133,7 +133,7 @@ def roll_func(x, axis=None): ts_pre.values, ts_pre_ano.values) # The previous step might create negative values (unlikely). Clip them - ts_pre.values = ts_pre.values.clip(0) + ts_pre.values = utils.clip_min(ts_pre.values, 0) assert np.all(np.isfinite(ts_pre.values)) assert np.all(np.isfinite(ts_tmp.values)) diff --git a/oggm/core/gis.py b/oggm/core/gis.py index b1dc34305..4e95fec66 100644 --- a/oggm/core/gis.py +++ b/oggm/core/gis.py @@ -52,7 +52,7 @@ pass # Locals -from oggm import entity_task +from oggm import entity_task, utils import oggm.cfg as cfg from oggm.exceptions import (InvalidParamsError, InvalidGeometryError, InvalidDEMError, GeometryError) @@ -356,7 +356,7 @@ def define_glacier_region(gdir, entity=None): .format(dxmethod)) # Additional trick for varying dx if dxmethod in ['linear', 'square']: - dx = np.clip(dx, cfg.PARAMS['d2'], cfg.PARAMS['dmax']) + dx = utils.clip_scalar(dx, cfg.PARAMS['d2'], cfg.PARAMS['dmax']) log.debug('(%s) area %.2f km, dx=%.1f', gdir.rgi_id, area, dx) @@ -573,7 +573,7 @@ def glacier_masks(gdir): dem_dr.close() # Clip topography to 0 m a.s.l. - dem = dem.clip(0) + dem = utils.clip_min(dem, 0) # Smooth DEM? if cfg.PARAMS['smooth_window'] > 0.: @@ -774,7 +774,7 @@ def simple_glacier_masks(gdir): dem_dr.close() # Clip topography to 0 m a.s.l. - dem = dem.clip(0) + dem = utils.clip_min(dem, 0) # Smooth DEM? if cfg.PARAMS['smooth_window'] > 0.: @@ -958,8 +958,9 @@ def gridded_attributes(gdir): glen_n = cfg.PARAMS['glen_n'] sy, sx = np.gradient(topo_smoothed, dx, dx) slope = np.arctan(np.sqrt(sy**2 + sx**2)) - slope_factor = np.clip(slope, np.deg2rad(cfg.PARAMS['min_slope']*4), - np.pi/2) + slope_factor = utils.clip_scalar(slope, + np.deg2rad(cfg.PARAMS['min_slope']*4), + np.pi/2) slope_factor = 1 / slope_factor**(glen_n / (glen_n+2)) aspect = np.arctan2(-sx, sy) @@ -1306,7 +1307,7 @@ def merged_glacier_masks(gdir, geometry): dem_dr.close() # Clip topography to 0 m a.s.l. - dem = dem.clip(0) + dem = utils.clip_min(dem, 0) # Interpolate shape to a regular path glacier_poly_hr = tolist(geometry) diff --git a/oggm/core/inversion.py b/oggm/core/inversion.py index 007ac318b..e2471b4e6 100644 --- a/oggm/core/inversion.py +++ b/oggm/core/inversion.py @@ -86,7 +86,7 @@ def prepare_for_inversion(gdir, add_debug_var=False, # Clip flux to 0 if np.any(flux < -0.1): log.warning('(%s) has negative flux somewhere', gdir.rgi_id) - flux = flux.clip(0) + flux = utils.clip_min(flux, 0) if fl.flows_to is None and gdir.inversion_calving_rate == 0: if not np.allclose(flux[-1], 0., atol=0.1): @@ -535,7 +535,7 @@ def distribute_thickness_per_altitude(gdir, add_slope=True, thick = np.where(glacier_mask, thick, 0.) # Re-mask - thick = thick.clip(0) + thick = utils.clip_min(thick, 0) thick[glacier_mask == 0] = np.NaN assert np.all(np.isfinite(thick[glacier_mask == 1])) @@ -625,7 +625,7 @@ def distribute_thickness_interp(gdir, add_slope=True, smooth_radius=None, points = np.array((np.ravel(yy[pok]), np.ravel(xx[pok]))).T inter = np.array((np.ravel(yy[pnan]), np.ravel(xx[pnan]))).T thick[pnan] = griddata(points, np.ravel(thick[pok]), inter, method='cubic') - thick = thick.clip(0) + thick = utils.clip_min(thick, 0) # Slope thick *= slope_factor @@ -704,7 +704,7 @@ def calving_flux_from_depth(gdir, k=None, water_depth=None, thick=None, fl = gdir.read_pickle('inversion_flowlines')[-1] # Altitude at the terminus and frontal width - t_altitude = np.clip(fl.surface_h[-1], 0, None) + t_altitude = utils.clip_min(fl.surface_h[-1], 0) width = fl.widths[-1] * gdir.grid.dx # Calving formula @@ -723,7 +723,7 @@ def calving_flux_from_depth(gdir, k=None, water_depth=None, thick=None, # Recompute free board before returning t_altitude = thick - water_depth - return {'flux': np.clip(flux, 0, None), + return {'flux': utils.clip_min(flux, 0), 'width': width, 'thick': thick, 'water_depth': water_depth, @@ -769,7 +769,7 @@ def find_inversion_calving_loop(gdir, initial_water_depth=None, max_ite=30, # Input if initial_water_depth is None: fl = gdir.read_pickle('inversion_flowlines')[-1] - initial_water_depth = np.clip(fl.surface_h[-1] / 3, 10, None) + initial_water_depth = utils.clip_min(fl.surface_h[-1] / 3, 10) rho = cfg.PARAMS['ice_density'] diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index 6e0dfd16a..6d6a45443 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -10,7 +10,7 @@ from oggm.cfg import SEC_IN_YEAR, SEC_IN_MONTH from oggm.utils import (SuperclassMeta, lazy_property, floatyear_to_date, date_to_floatyear, monthly_timeseries, ncDataset, - tolist) + tolist, clip_min, clip_max) class MassBalanceModel(object, metaclass=SuperclassMeta): @@ -210,7 +210,7 @@ def temp_bias(self, value): def get_monthly_mb(self, heights, year=None, fl_id=None): mb = (np.asarray(heights) - self.ela_h) * self.grad if self.max_mb is not None: - mb = mb.clip(None, self.max_mb) + mb = clip_max(mb, self.max_mb) return mb / SEC_IN_YEAR / self.rho def get_annual_mb(self, heights, year=None, fl_id=None): @@ -373,7 +373,7 @@ def get_monthly_climate(self, heights, year=None): npix = len(heights) temp = np.ones(npix) * itemp + igrad * (heights - self.ref_hgt) tempformelt = temp - self.t_melt - tempformelt[:] = np.clip(tempformelt, 0., None) + tempformelt[:] = clip_min(tempformelt, 0.) # Compute solid precipitation from total precipitation prcp = np.ones(npix) * iprcp @@ -408,13 +408,12 @@ def _get_2d_annual_climate(self, heights, year): self.ref_hgt) temp2d = np.atleast_2d(itemp).repeat(npix, 0) + grad_temp temp2dformelt = temp2d - self.t_melt - temp2dformelt[:] = np.clip(temp2dformelt, 0, None) + temp2dformelt[:] = clip_min(temp2dformelt, 0) # Compute solid precipitation from total precipitation prcp = np.atleast_2d(iprcp).repeat(npix, 0) fac = 1 - (temp2d - self.t_solid) / (self.t_liq - self.t_solid) - fac = np.clip(fac, 0, 1) - prcpsol = prcp * fac + prcpsol = prcp * np.clip(fac, 0, 1) return temp2d, temp2dformelt, prcp, prcpsol diff --git a/oggm/core/sia2d.py b/oggm/core/sia2d.py index 1c8dff72d..a2d54829d 100644 --- a/oggm/core/sia2d.py +++ b/oggm/core/sia2d.py @@ -421,10 +421,12 @@ def step(self, dt): div_q, dt_cfl = self.diffusion_upstream_2d() - dt_use = np.clip(np.min([dt_cfl, dt]), 0, self.max_dt) + dt_use = utils.clip_scalar(np.min([dt_cfl, dt]), 0, self.max_dt) - self.ice_thick = (self.surface_h + (self.get_mb() + div_q) * dt_use - - self.bed_topo).clip(0) + self.ice_thick = utils.clip_min(self.surface_h + + (self.get_mb() + div_q) * dt_use - + self.bed_topo, + 0) # Next step self.t += dt_use diff --git a/oggm/tests/test_utils.py b/oggm/tests/test_utils.py index 6450581bd..c0b88d55a 100644 --- a/oggm/tests/test_utils.py +++ b/oggm/tests/test_utils.py @@ -251,6 +251,57 @@ def test_adhikari_shape_factors(self): np.testing.assert_equal(shape_factor_adhikari(w, h, is_rect*0.), factors_parabolic) + def test_clips(self): + + a = np.arange(50) - 5 + assert_array_equal(np.clip(a, 0, 10), utils.clip_scalar(a, 0, 10)) + assert_array_equal(np.clip(a, None, 10), utils.clip_max(a, 10)) + assert_array_equal(np.clip(a, 0, None), utils.clip_min(a, 0)) + + for a in [-1, 2, 12]: + assert_array_equal(np.clip(a, 0, 10), utils.clip_scalar(a, 0, 10)) + assert_array_equal(np.clip(a, None, 10), utils.clip_max(a, 10)) + assert_array_equal(np.clip(a, 0, None), utils.clip_min(a, 0)) + + def test_clips_performance(self): + import timeit + + n = int(1e5) + + # Scalar + s1 = 'import numpy as np' + s2 = 'from oggm import utils' + + t1 = timeit.timeit('np.clip(1, 0, 10)', number=n, setup=s1) + t2 = timeit.timeit('utils.clip_scalar(1, 0, 10)', number=n, setup=s2) + assert t2 < t1 + + t1 = timeit.timeit('np.clip(12, None, 10)', number=n, setup=s1) + t2 = timeit.timeit('utils.clip_max(12, 10)', number=n, setup=s2) + assert t2 < t1 + + t1 = timeit.timeit('np.clip(12, 15, None)', number=n, setup=s1) + t2 = timeit.timeit('utils.clip_min(12, 15)', number=n, setup=s2) + assert t2 < t1 + + # Array + s1 = 'import numpy as np; a = np.arange(50) - 5' + s2 = 'from oggm import utils; import numpy as np; a = np.arange(50)-5' + + t1 = timeit.timeit('np.clip(a, 0, 10)', number=n, setup=s1) + t2 = timeit.timeit('utils.clip_scalar(a, 0, 10)', number=n, setup=s2) + # This usually fails as advertised by numpy + # (although with np 1.17 I'm not sure) + # assert t2 < t1 + + t1 = timeit.timeit('np.clip(a, None, 10)', number=n, setup=s1) + t2 = timeit.timeit('utils.clip_max(a, 10)', number=n, setup=s2) + assert t2 < t1 + + t1 = timeit.timeit('np.clip(a, 15, None)', number=n, setup=s1) + t2 = timeit.timeit('utils.clip_min(a, 15)', number=n, setup=s2) + assert t2 < t1 + class TestInitialize(unittest.TestCase): diff --git a/oggm/utils/_funcs.py b/oggm/utils/_funcs.py index 099def83e..19fb41d8b 100644 --- a/oggm/utils/_funcs.py +++ b/oggm/utils/_funcs.py @@ -333,6 +333,21 @@ def corrcoef(ref, data): return np.corrcoef(ref, data)[0, 1] +def clip_scalar(value, vmin, vmax): + """A faster numpy.clip ON SCALARS ONLY. + + See https://github.com/numpy/numpy/issues/14281 + """ + return np.core.umath.maximum(np.core.umath.minimum(value, vmax), vmin) + + +# A faster numpy.clip when only one value is clipped (here: min). +clip_min = np.core.umath.maximum + +# A faster numpy.clip when only one value is clipped (here: max). +clip_max = np.core.umath.minimum + + def nicenumber(number, binsize, lower=False): """Returns the next higher or lower "nice number", given by binsize. From f3cd0103d275408282f9011cdb96af890f537a24 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 26 Sep 2019 12:30:52 +0100 Subject: [PATCH 2/2] Further optims --- oggm/core/gis.py | 6 +++--- oggm/core/inversion.py | 6 +++--- oggm/core/massbalance.py | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/oggm/core/gis.py b/oggm/core/gis.py index 4e95fec66..8921a7c76 100644 --- a/oggm/core/gis.py +++ b/oggm/core/gis.py @@ -573,7 +573,7 @@ def glacier_masks(gdir): dem_dr.close() # Clip topography to 0 m a.s.l. - dem = utils.clip_min(dem, 0) + utils.clip_min(dem, 0, out=dem) # Smooth DEM? if cfg.PARAMS['smooth_window'] > 0.: @@ -774,7 +774,7 @@ def simple_glacier_masks(gdir): dem_dr.close() # Clip topography to 0 m a.s.l. - dem = utils.clip_min(dem, 0) + utils.clip_min(dem, 0, out=dem) # Smooth DEM? if cfg.PARAMS['smooth_window'] > 0.: @@ -1307,7 +1307,7 @@ def merged_glacier_masks(gdir, geometry): dem_dr.close() # Clip topography to 0 m a.s.l. - dem = utils.clip_min(dem, 0) + utils.clip_min(dem, 0, out=dem) # Interpolate shape to a regular path glacier_poly_hr = tolist(geometry) diff --git a/oggm/core/inversion.py b/oggm/core/inversion.py index e2471b4e6..5790e8b14 100644 --- a/oggm/core/inversion.py +++ b/oggm/core/inversion.py @@ -86,7 +86,7 @@ def prepare_for_inversion(gdir, add_debug_var=False, # Clip flux to 0 if np.any(flux < -0.1): log.warning('(%s) has negative flux somewhere', gdir.rgi_id) - flux = utils.clip_min(flux, 0) + utils.clip_min(flux, 0, out=flux) if fl.flows_to is None and gdir.inversion_calving_rate == 0: if not np.allclose(flux[-1], 0., atol=0.1): @@ -535,7 +535,7 @@ def distribute_thickness_per_altitude(gdir, add_slope=True, thick = np.where(glacier_mask, thick, 0.) # Re-mask - thick = utils.clip_min(thick, 0) + utils.clip_min(thick, 0, out=thick) thick[glacier_mask == 0] = np.NaN assert np.all(np.isfinite(thick[glacier_mask == 1])) @@ -625,7 +625,7 @@ def distribute_thickness_interp(gdir, add_slope=True, smooth_radius=None, points = np.array((np.ravel(yy[pok]), np.ravel(xx[pok]))).T inter = np.array((np.ravel(yy[pnan]), np.ravel(xx[pnan]))).T thick[pnan] = griddata(points, np.ravel(thick[pok]), inter, method='cubic') - thick = utils.clip_min(thick, 0) + utils.clip_min(thick, 0, out=thick) # Slope thick *= slope_factor diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index 6d6a45443..c5355a624 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -210,7 +210,7 @@ def temp_bias(self, value): def get_monthly_mb(self, heights, year=None, fl_id=None): mb = (np.asarray(heights) - self.ela_h) * self.grad if self.max_mb is not None: - mb = clip_max(mb, self.max_mb) + clip_max(mb, self.max_mb, out=mb) return mb / SEC_IN_YEAR / self.rho def get_annual_mb(self, heights, year=None, fl_id=None): @@ -373,7 +373,7 @@ def get_monthly_climate(self, heights, year=None): npix = len(heights) temp = np.ones(npix) * itemp + igrad * (heights - self.ref_hgt) tempformelt = temp - self.t_melt - tempformelt[:] = clip_min(tempformelt, 0.) + clip_min(tempformelt, 0, out=tempformelt) # Compute solid precipitation from total precipitation prcp = np.ones(npix) * iprcp @@ -408,7 +408,7 @@ def _get_2d_annual_climate(self, heights, year): self.ref_hgt) temp2d = np.atleast_2d(itemp).repeat(npix, 0) + grad_temp temp2dformelt = temp2d - self.t_melt - temp2dformelt[:] = clip_min(temp2dformelt, 0) + clip_min(temp2dformelt, 0, out=temp2dformelt) # Compute solid precipitation from total precipitation prcp = np.atleast_2d(iprcp).repeat(npix, 0)