Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# Allow parabolic beds to grow
mb = dt * mb * np.where((mb > 0.) & (widths == 0), 10., widths)
# Update section with ice flow and mass balance
new_section = (fl.section + (flx_stag[0:-1] - flx_stag[1:])*dt/dx +
trib_flux*dt/dx + mb)
# Keep positive values only and store
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]] += \
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 += utils.clip_min(flx_stag[-2], 0)*dt
# Next step
self.t += dt
return dt
assert len(bed) == idl.nx
pvalid = np.sum(np.isfinite(bed)) / len(bed) * 100
log.debug('(%s) percentage of valid parabolas in downstream: %d',
gdir.rgi_id, int(pvalid))
# Scale for dx (we worked in grid coords but need meters)
bed = bed / gdir.grid.dx**2
# interpolation, filling the gaps
default = cfg.PARAMS['default_parabolic_bedshape']
bed_int = interp_nans(bed, default=default)
# 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 = utils.clip_min(bed_int, cfg.PARAMS['downstream_min_shape'])
# Smoothing
bed_ma = pd.Series(bed_int)
bed_ma = bed_ma.rolling(window=5, center=True, min_periods=1).mean()
return bed_ma.values
if thick is None:
thick = cl['thick'][-1]
if water_depth is None:
water_depth = thick - t_altitude
elif not fixed_water_depth:
# Correct thickness with prescribed water depth
# If fixed_water_depth=True then we forget about t_altitude
thick = water_depth + t_altitude
flux = k * thick * water_depth * width / 1e9
if fixed_water_depth:
# Recompute free board before returning
t_altitude = thick - water_depth
return {'flux': utils.clip_min(flux, 0),
'width': width,
'thick': thick,
'water_depth': water_depth,
'free_board': t_altitude}
# Default gradient?
if igrad is None:
igrad = itemp * 0 + default_grad
# Correct precipitation
iprcp *= prcp_fac
# For each height pixel:
# Compute temp and tempformelt (temperature above melting threshold)
npix = len(heights)
grad_temp = np.atleast_2d(igrad).repeat(npix, 0)
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 = 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)
fac = utils.clip_array(fac, 0, 1)
prcpsol = prcpsol * fac
return time, temp2dformelt, prcpsol
dis_from_border /= np.mean(dis_from_border[glacier_mask == 1])
thick *= dis_from_border
# Slope
thick *= slope_factor
# Smooth
dx = gdir.grid.dx
if smooth_radius != 0:
if smooth_radius is None:
smooth_radius = np.rint(cfg.PARAMS['smooth_window'] / dx)
thick = gaussian_blur(thick, np.int(smooth_radius))
thick = np.where(glacier_mask, thick, 0.)
# Re-mask
utils.clip_min(thick, 0, out=thick)
thick[glacier_mask == 0] = np.NaN
assert np.all(np.isfinite(thick[glacier_mask == 1]))
# Conserve volume
tmp_vol = np.nansum(thick * dx**2)
thick *= init_vol / tmp_vol
# write
with utils.ncDataset(grids_file, 'a') as nc:
vn = 'distributed_thickness' + varname_suffix
if vn in nc.variables:
v = nc.variables[vn]
else:
v = nc.createVariable(vn, 'f4', ('y', 'x', ), zlib=True)
v.units = '-'
v.long_name = 'Distributed ice thickness'
# For these we had an additional grid point
if is_trib or self.is_tidewater:
flx_stag = flx_stag[:-1]
# Mass balance
widths = fl.widths_m
mb = self.get_mb(fl.surface_h, self.yr, fl_id=fl_id)
# Allow parabolic beds to grow
mb = dt * mb * np.where((mb > 0.) & (widths == 0), 10., widths)
# Update section with ice flow and mass balance
new_section = (fl.section + (flx_stag[0:-1] - flx_stag[1:])*dt/dx +
trib_flux*dt/dx + mb)
# Keep positive values only and store
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]] += \
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 += utils.clip_min(flx_stag[-2], 0)*dt
# Next step
self.t += dt
return dt
# Widths
widths = fl.widths * gdir.grid.dx
# Heights
hgt = fl.surface_h
angle = -np.gradient(hgt, dx) # beware the minus sign
# Flux needs to be in [m3 s-1] (*ice* velocity * surface)
# fl.flux is given in kg m-2 yr-1, rho in kg m-3, so this should be it:
rho = cfg.PARAMS['ice_density']
flux = fl.flux * (gdir.grid.dx**2) / cfg.SEC_IN_YEAR / rho
# Clip flux to 0
if np.any(flux < -0.1):
log.warning('(%s) has negative flux somewhere', gdir.rgi_id)
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):
# TODO: this test doesn't seem meaningful here
msg = ('({}) flux at terminus should be zero, but is: '
'{.4f} m3 ice s-1'.format(gdir.rgi_id, flux[-1]))
raise RuntimeError(msg)
flux[-1] = 0.
# Shape
is_rectangular = fl.is_rectangular
if not invert_with_rectangular:
is_rectangular[:] = False
if invert_all_rectangular:
is_rectangular[:] = True
log.info('(%s) local mu* computation for t*=%d', gdir.rgi_id, tstar)
# Get the corresponding mu
years, temp_yr, prcp_yr = mb_yearly_climate_on_glacier(gdir, year_range=yr)
assert len(years) == (2 * mu_hp + 1)
# mustar is taking calving into account (units of specific MB)
mustar = (np.mean(prcp_yr) - cmb) / np.mean(temp_yr)
if not np.isfinite(mustar):
raise MassBalanceCalibrationError('{} has a non finite '
'mu'.format(gdir.rgi_id))
# Clip it?
if cfg.PARAMS['clip_mu_star']:
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']):
raise MassBalanceCalibrationError('mu* out of specified bounds: '
'{:.2f}'.format(mustar))
# Scalars in a small dict for later
df = dict()
df['rgi_id'] = gdir.rgi_id
df['t_star'] = int(tstar)
df['bias'] = bias
df['mu_star_glacierwide'] = mustar
gdir.write_json(df, 'local_mustar')
# Param
nmin = int(cfg.PARAMS['min_n_per_bin'])
smooth_ws = int(cfg.PARAMS['smooth_widths_window_size'])
# Per flowline (important so that later, the indices can be moved)
catchment_heights = []
for ci in catchment_indices:
_t = topo[tuple(ci.T)][:]
catchment_heights.append(list(_t[np.isfinite(_t)]))
# Loop over lines in a reverse order
for fl, catch_h in zip(fls, catchment_heights):
# Interpolate widths
widths = utils.interp_nans(fl.widths)
widths = utils.clip_min(widths, 0.1)
# Get topo per catchment and per flowline point
fhgt = fl.surface_h
# Max and mins for the histogram
maxh = np.max(fhgt)
minh = np.min(fhgt)
# Sometimes, the centerline does not reach as high as each pix on the
# glacier. (e.g. RGI40-11.00006)
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 = utils.clip_min(catch_h, minh)