Use faster clips everywhere (#873)

* Use faster clips everywhere

* Further optims

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']

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']):

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']

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))

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)
+    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 = dem.clip(0)
+    utils.clip_min(dem, 0, out=dem)
 
     # 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)
+    utils.clip_min(dem, 0, out=dem)
 
     # Interpolate shape to a regular path
     glacier_poly_hr = tolist(geometry)

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)
+        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 = thick.clip(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 = thick.clip(0)
+    utils.clip_min(thick, 0, out=thick)
 
     # 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']
 

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)
+            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[:] = np.clip(tempformelt, 0., None)
+        clip_min(tempformelt, 0, out=tempformelt)
 
         # 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)
+        clip_min(temp2dformelt, 0, out=temp2dformelt)
 
         # 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
 

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