Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
diag_n_bad_slopes = 0
diag_n_pix = 0
for ic, cl in enumerate(cls):
points = line_interpol(cl.line, dx)
# For tributaries, remove the tail
if ic < (len(cls)-1):
points = points[0:-lid]
new_line = shpg.LineString(points)
# Interpolate heights
xx, yy = new_line.xy
hgts = interpolator((yy, xx))
if len(hgts) < 5:
raise GeometryError('This centerline is too short')
# If smoothing, this is the moment
hgts = gaussian_filter1d(hgts, sw)
# Check for min slope issues and correct if needed
if do_filter:
# Correct only where glacier
hgts = _filter_small_slopes(hgts, dx*gdir.grid.dx)
isfin = np.isfinite(hgts)
if not np.any(isfin):
raise GeometryError('This centerline has no positive slopes')
diag_n_bad_slopes += np.sum(~isfin)
diag_n_pix += len(isfin)
perc_bad = np.sum(~isfin) / len(isfin)
if perc_bad > 0.8:
log.warning('({}) more than {:.0%} of the flowline is cropped '
# Catchment polygon
mask[:] = 0
mask[tuple(ci.T)] = 1
poly, poly_no = _mask_to_polygon(mask, gdir=gdir)
# First guess widths
for i, (normal, pcoord) in enumerate(zip(fl.normals, fl.line.coords)):
width, wline = _point_width(normal, pcoord, fl, poly, poly_no)
widths[i] = width
wlines.append(wline)
valid = np.where(np.isfinite(widths))
if len(valid[0]) == 0:
errmsg = '({}) first guess widths went wrong.'.format(gdir.rgi_id)
raise GeometryError(errmsg)
# Ok now the entire centerline is computed.
# I take all these widths for geometrically valid, and see if they
# intersect with our buffered catchment/glacier intersections
is_rectangular = []
for wg in wlines:
is_rectangular.append(np.any(gdfi.intersects(wg)))
is_rectangular = _filter_grouplen(is_rectangular, minsize=5)
# we filter the lines which have a large altitude range
fil_widths = _filter_for_altitude_range(widths, wlines, topo)
# Filter +- widths at junction points
for fid in fl.inflow_indices:
i0 = int(utils.clip_scalar(fid-jpix, jpix/2, n-jpix/2))
i1 = int(utils.clip_scalar(fid+jpix+1, jpix/2, n-jpix/2))
# And rejoin the cutted tails
olines = _join_lines(olines, oheads)
# Adds the line level
for cl in olines:
cl.order = line_order(cl)
# And sort them per order !!! several downstream tasks rely on this
cls = []
for i in np.argsort([cl.order for cl in olines]):
cls.append(olines[i])
# Final check
if len(cls) == 0:
raise GeometryError('({}) no valid centerline could be '
'found!'.format(gdir.rgi_id))
# Write the data
gdir.write_pickle(cls, 'centerlines')
if is_first_call:
# For diagnostics of filtered centerlines
gdir.add_to_diagnostics('n_orig_centerlines', len(cls))
def glacier_masks(gdir):
"""Makes a gridded mask of the glacier outlines that can be used by OGGM.
For a more robust solution (not OGGM compatible) see simple_glacier_masks.
Parameters
----------
gdir : :py:class:`oggm.GlacierDirectory`
where to write the data
"""
# In case nominal, just raise
if gdir.is_nominal:
raise GeometryError('{} is a nominal glacier.'.format(gdir.rgi_id))
if not os.path.exists(gdir.get_filepath('gridded_data')):
# In a possible future, we might actually want to raise a
# deprecation warning here
process_dem(gdir)
# Geometries
geometry = gdir.read_shapefile('outlines').geometry[0]
# Interpolate shape to a regular path
glacier_poly_hr = _interp_polygon(geometry, gdir.grid.dx)
# Transform geometry into grid coordinates
# It has to be in pix center coordinates because of how skimage works
def proj(x, y):
grid = gdir.grid.center_grid
xc.extend(dwl[:, 0])
yc.extend(dwl[:, 1])
altrange = topo[yc, xc]
if len(np.where(np.isfinite(altrange))[0]) != 0:
if (np.nanmax(altrange) - np.nanmin(altrange)) > alt_range_th:
out_width[i] = np.NaN
valid = np.where(np.isfinite(out_width))
if len(valid[0]) > 0:
break
else:
alt_range_th += 20
log.warning('Set altitude threshold to {}'.format(alt_range_th))
if alt_range_th > 2000:
raise GeometryError('Problem by altitude filter.')
return out_width
# to the last longest line
diff = l.difference(toremove)
if diff.type is 'MultiLineString':
# Remove the lines that have no head
diff = list(diff)
for il in diff:
hashead = False
for h in heads:
if il.intersects(h):
hashead = True
diff = il
break
if hashead:
break
else:
raise GeometryError('Head not found')
# keep this head line only if it's long enough
if diff.length > r:
# Fun fact. The heads can be cut by the buffer too
diff = shpg.LineString(l.coords[0:2] + diff.coords[2:])
tokeep.append(diff)
ilines = tokeep
# it could happen that we're done at this point
if len(ilines) == 0:
break
# Otherwise keep the longest one and continue
lengths = np.array([])
for l in ilines:
lengths = np.append(lengths, l.length)
ll = ilines[np.argmax(lengths)]
if (ss / sr) > 0.2:
log.warning('(%s) this blob was unusually large', rid)
mask[:] = 0
mask[np.where(regions == (am+1))] = 1
nlist = measure.find_contours(mask, 0.5)
# First is the exterior, the rest are nunataks
e_line = shpg.LinearRing(nlist[0][:, ::-1])
i_lines = [shpg.LinearRing(ipoly[:, ::-1]) for ipoly in nlist[1:]]
poly = shpg.Polygon(e_line, i_lines).buffer(0)
if not poly.is_valid:
raise GeometryError('Mask to polygon conversion error.')
poly_no = shpg.Polygon(e_line).buffer(0)
if not poly_no.is_valid:
raise GeometryError('Mask to polygon conversion error.')
return poly, poly_no