How to use the oggm.exceptions.GeometryError function in oggm

To help you get started, we’ve selected a few oggm examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github OGGM / oggm / oggm / core / centerlines.py View on Github external
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 '
github OGGM / oggm / oggm / core / centerlines.py View on Github external
# 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))
github OGGM / oggm / oggm / core / centerlines.py View on Github external
# 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))
github OGGM / oggm / oggm / core / gis.py View on Github external
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
github OGGM / oggm / oggm / core / centerlines.py View on Github external
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
github OGGM / oggm / oggm / core / centerlines.py View on Github external
# 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)]
github OGGM / oggm / oggm / core / centerlines.py View on Github external
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