Add diagnostic tools which make it easier to detect dubious glacier o…

…utlines (#445)

* Add some RGI tools

* Add simple_glacier_mask task

* What's new and test failure

* Add flowline filters to diagnostics

* Color docs

* Add n_orig_centerlines to diagnostics

docs/_templates/layout.html

@@ -4,7 +4,7 @@
   <style>
     /* Sidebar header (and topbar for mobile) */
     .wy-side-nav-search, .wy-nav-top {
-      background: #9eceed;
+      background: #7096b0;
     }
   </style>
 {% endblock %}

docs/whats-new.rst

@@ -60,6 +60,10 @@ Enhancements
   along with the original data citation (:pull:`441`). We encourage
   to always cite the original data providers.
   By `Fabien Maussion <https://github.com/fmaussion>`_.
+- Added some diagnostic tools which make it easier to detect dubious glacier
+  outlines or DEM errors (:pull:`445`). This will be used to report to the
+  RGI authors.
+  By `Fabien Maussion <https://github.com/fmaussion>`_.
 
 
 Bug fixes

oggm/cfg.py

@@ -127,6 +127,9 @@ def __setitem__(self, key, value):
 _doc = 'A ``salem.Grid`` handling the georeferencing of the local grid.'
 BASENAMES['glacier_grid'] = ('glacier_grid.json', _doc)
 
+_doc = 'A dictionary containing runtime diagnostics useful for debugging.'
+BASENAMES['diagnostics'] = ('diagnostics.json', _doc)
+
 _doc = 'A netcdf file containing several gridded data variables such as ' \
        'topography, the glacier masks and more (see the netCDF file metadata).'
 BASENAMES['gridded_data'] = ('gridded_data.nc', _doc)

oggm/core/centerlines.py

@@ -806,7 +806,11 @@ def compute_centerlines(gdir, heads=None):
     ext_yx = tuple(reversed(poly_pix.exterior.xy))
     zoutline = topo[y[:-1], x[:-1]]  # last point is first point
 
+    # For diagnostics
+    is_first_call = False
     if heads is None:
+        # This is the default for when no filter is yet applied
+        is_first_call = True
         heads = _get_centerlines_heads(gdir, ext_yx, zoutline, single_fl,
                                        glacier_mask, topo, geom, poly_pix)
 
@@ -857,6 +861,10 @@ def compute_centerlines(gdir, heads=None):
     # 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))
+
     # Netcdf
     with netCDF4.Dataset(grids_file, 'a') as nc:
         if 'cost_grid' in nc.variables:
@@ -1551,7 +1559,8 @@ def initialize_flowlines(gdir):
 
     # Smooth window
     sw = cfg.PARAMS['flowline_height_smooth']
-
+    diag_n_bad_slopes = 0
+    diag_n_pix = 0
     for ic, cl in enumerate(cls):
         points = line_interpol(cl.line, dx)
 
@@ -1574,7 +1583,10 @@ def initialize_flowlines(gdir):
             # Correct only where glacier
             hgts = _filter_small_slopes(hgts, dx*gdir.grid.dx)
             isfin = np.isfinite(hgts)
-            assert np.any(isfin)
+            if not np.any(isfin):
+                raise RuntimeError('This line 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 '
@@ -1602,6 +1614,9 @@ def initialize_flowlines(gdir):
 
     # Write the data
     gdir.write_pickle(fls, 'inversion_flowlines')
+    if do_filter:
+        out = diag_n_bad_slopes/diag_n_pix
+        gdir.add_to_diagnostics('perc_invalid_flowline', out)
 
 
 @entity_task(log, writes=['inversion_flowlines'])

oggm/core/gis.py

@@ -189,10 +189,10 @@ def _polygon_to_pix(polygon):
                     if tmp.is_valid:
                         break
                 if b == 0.99:
-                    raise RuntimeError('This glacier geometry is crazy.')
+                    raise RuntimeError('This glacier geometry is not valid.')
 
     if not tmp.is_valid:
-        raise RuntimeError('This glacier geometry is crazy.')
+        raise RuntimeError('This glacier geometry is not valid.')
 
     return tmp
 
@@ -469,7 +469,7 @@ def proj(x, y):
     # fix-invalid-polygon-python-shapely
     glacier_poly_hr = glacier_poly_hr.buffer(0)
     if not glacier_poly_hr.is_valid:
-        raise RuntimeError('This glacier geometry is crazy.')
+        raise RuntimeError('This glacier geometry is not valid.')
 
     # Rounded nearest pix
     glacier_poly_pix = _polygon_to_pix(glacier_poly_hr)
@@ -548,3 +548,167 @@ def proj(x, y):
     geometries['polygon_pix'] = glacier_poly_pix
     geometries['polygon_area'] = geometry.area
     gdir.write_pickle(geometries, 'geometries')
+
+
+@entity_task(log, writes=['gridded_data'])
+def simple_glacier_masks(gdir):
+    """Compute glacier masks based on much simpler rules than the default.
+
+    This is therefore more robust
+
+    Parameters
+    ----------
+    gdir : :py:class:`oggm.GlacierDirectory`
+        where to write the data
+    """
+
+    # open srtm tif-file:
+    dem_dr = rasterio.open(gdir.get_filepath('dem'), 'r', driver='GTiff')
+    dem = dem_dr.read(1).astype(rasterio.float32)
+
+    # Grid
+    nx = dem_dr.width
+    ny = dem_dr.height
+    assert nx == gdir.grid.nx
+    assert ny == gdir.grid.ny
+
+    # Correct the DEM (ASTER...)
+    # Currently we just do a linear interp -- ASTER is totally shit anyway
+    min_z = -999.
+    isfinite = np.isfinite(dem)
+    if (np.min(dem) <= min_z) or np.any(~isfinite):
+        xx, yy = gdir.grid.ij_coordinates
+        pnan = np.nonzero((dem <= min_z) | (~isfinite))
+        pok = np.nonzero((dem > min_z) | isfinite)
+        points = np.array((np.ravel(yy[pok]), np.ravel(xx[pok]))).T
+        inter = np.array((np.ravel(yy[pnan]), np.ravel(xx[pnan]))).T
+        dem[pnan] = griddata(points, np.ravel(dem[pok]), inter)
+        log.warning(gdir.rgi_id + ': DEM needed interpolation.')
+
+    isfinite = np.isfinite(dem)
+    if not np.all(isfinite):
+        # see how many percent of the dem
+        if np.sum(~isfinite) > (0.2 * nx * ny):
+            raise RuntimeError('({}) too many NaNs in DEM'.format(gdir.rgi_id))
+        log.warning('({}) DEM needed zeros somewhere.'.format(gdir.rgi_id))
+        dem[isfinite] = 0
+
+    if np.min(dem) == np.max(dem):
+        raise RuntimeError('({}) min equal max in the DEM.'
+                           .format(gdir.rgi_id))
+
+    # Proj
+    if LooseVersion(rasterio.__version__) >= LooseVersion('1.0'):
+        transf = dem_dr.transform
+    else:
+        transf = dem_dr.affine
+    x0 = transf[2]  # UL corner
+    y0 = transf[5]  # UL corner
+    dx = transf[0]
+    dy = transf[4]  # Negative
+    assert dx == -dy
+    assert dx == gdir.grid.dx
+    assert y0 == gdir.grid.corner_grid.y0
+    assert x0 == gdir.grid.corner_grid.x0
+    dem_dr.close()
+
+    # Clip topography to 0 m a.s.l.
+    dem = dem.clip(0)
+
+    # Smooth DEM?
+    if cfg.PARAMS['smooth_window'] > 0.:
+        gsize = np.rint(cfg.PARAMS['smooth_window'] / dx)
+        smoothed_dem = gaussian_blur(dem, np.int(gsize))
+    else:
+        smoothed_dem = dem.copy()
+
+    if not np.all(np.isfinite(smoothed_dem)):
+        raise RuntimeError('({}) NaN in smoothed DEM'.format(gdir.rgi_id))
+
+    # Geometries
+    outlines_file = gdir.get_filepath('outlines')
+    geometry = gpd.GeoDataFrame.from_file(outlines_file).geometry[0]
+
+    # 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
+        return grid.transform(x, y, crs=grid.proj)
+    geometry = shapely.ops.transform(proj, geometry)
+
+    # simple trick to correct invalid polys:
+    # http://stackoverflow.com/questions/20833344/
+    # fix-invalid-polygon-python-shapely
+    geometry = geometry.buffer(0)
+    if not geometry.is_valid:
+        raise RuntimeError('This glacier geometry is not valid.')
+
+    # Compute the glacier mask (currently: center pixels + touched)
+    nx, ny = gdir.grid.nx, gdir.grid.ny
+    glacier_mask = np.zeros((ny, nx), dtype=np.uint8)
+    glacier_ext = np.zeros((ny, nx), dtype=np.uint8)
+    x, y = geometry.exterior.xy
+    x, y = np.array(x), np.array(y)
+    glacier_mask[skdraw.polygon(y, x)] = 1
+    glacier_mask[skdraw.polygon_perimeter(y, x)] = 1
+    glacier_ext[skdraw.polygon_perimeter(y, x)] = 1
+    for gint in geometry.interiors:
+        x, y = gint.xy
+        x, y = np.array(x), np.array(y)
+        glacier_mask[skdraw.polygon(y, x)] = 0
+        glacier_mask[skdraw.polygon_perimeter(y, x)] = 0
+
+    # Because of the 0 values at nunataks boundaries, some "Ice Islands"
+    # can happen within nunataks (e.g.: RGI40-11.00062)
+    # See if we can filter them out easily
+    regions, nregions = label(glacier_mask, structure=label_struct)
+    if nregions > 1:
+        log.debug('(%s) we had to cut an island in the mask', gdir.rgi_id)
+        # Check the size of those
+        region_sizes = [np.sum(regions == r) for r in np.arange(1, nregions+1)]
+        am = np.argmax(region_sizes)
+        # Check not a strange glacier
+        sr = region_sizes.pop(am)
+        for ss in region_sizes:
+            assert (ss / sr) < 0.1
+        glacier_mask[:] = 0
+        glacier_mask[np.where(regions == (am+1))] = 1
+
+    # Last sanity check based on the masked dem
+    tmp_max = np.max(dem[np.where(glacier_mask == 1)])
+    tmp_min = np.min(dem[np.where(glacier_mask == 1)])
+    if tmp_max < (tmp_min + 1):
+        raise RuntimeError('({}) min equal max in the masked DEM.'
+                           .format(gdir.rgi_id))
+
+    # write out the grids in the netcdf file
+    nc = gdir.create_gridded_ncdf_file('gridded_data')
+
+    v = nc.createVariable('topo', 'f4', ('y', 'x', ), zlib=True)
+    v.units = 'm'
+    v.long_name = 'DEM topography'
+    v[:] = dem
+
+    v = nc.createVariable('topo_smoothed', 'f4', ('y', 'x', ), zlib=True)
+    v.units = 'm'
+    v.long_name = ('DEM topography smoothed' 
+                   ' with radius: {:.1} m'.format(cfg.PARAMS['smooth_window']))
+    v[:] = smoothed_dem
+
+    v = nc.createVariable('glacier_mask', 'i1', ('y', 'x', ), zlib=True)
+    v.units = '-'
+    v.long_name = 'Glacier mask'
+    v[:] = glacier_mask
+
+    v = nc.createVariable('glacier_ext', 'i1', ('y', 'x', ), zlib=True)
+    v.units = '-'
+    v.long_name = 'Glacier external boundaries'
+    v[:] = glacier_ext
+
+    # add some meta stats and close
+    nc.max_h_dem = np.max(dem)
+    nc.min_h_dem = np.min(dem)
+    dem_on_g = dem[np.where(glacier_mask)]
+    nc.max_h_glacier = np.max(dem_on_g)
+    nc.min_h_glacier = np.min(dem_on_g)
+    nc.close()

oggm/sandbox/rgi_tools/scripts.py

@@ -12,14 +12,14 @@
 from oggm.sandbox.rgi_tools.tools import (compute_intersects, prepare_divides,
                                           OUTDIR_INTERSECTS)
 
-rgi_version = '6'
+rgi_version = '61'
 
 
 def intersects_script():
     # Download RGI files
     cfg.initialize()
     rgi_dir = get_rgi_dir(version=rgi_version)
-    fp = '*_rgi' + rgi_version + '0_*.shp'
+    fp = '*_rgi' + rgi_version + '_*.shp'
     rgi_shps = list(glob(os.path.join(rgi_dir, "*", fp)))
     rgi_shps = [r for r in rgi_shps if 'Regions' not in r]
     with mp.Pool() as p:
@@ -72,4 +72,4 @@ def divides_script():
 
 
 if __name__ == "__main__":
-    merge_intersects()
+    intersects_script()

oggm/sandbox/rgi_tools/tools.py

@@ -15,7 +15,7 @@
 
 INDIR_DIVIDES = '/home/mowglie/disk/Data/OGGM_DATA/results_global_partitioning/altitude_filter/'
 
-OUTDIR_INTERSECTS = '/home/mowglie/tmp/RGI_V6_Intersects/'
+OUTDIR_INTERSECTS = '/home/mowglie/tmp/RGI_V61_Intersects/'
 OUTDIR_DIVIDES = '/home/mowglie/disk/Data/OGGM_DATA/RGI_V5_Modified/'
 
 

oggm/tasks.py

@@ -5,6 +5,7 @@
 # Entity tasks
 from oggm.core.gis import define_glacier_region
 from oggm.core.gis import glacier_masks
+from oggm.core.gis import simple_glacier_masks
 from oggm.core.centerlines import compute_centerlines
 from oggm.core.centerlines import compute_downstream_line
 from oggm.core.centerlines import compute_downstream_bedshape

oggm/tests/test_models.py

@@ -1756,6 +1756,7 @@ def test_inversion_noisy_sf_adhikari(self):
         cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = old_model_sf
         cfg.PARAMS['use_shape_factor_for_inversion'] = old_inversion_sf
 
+    @is_slow
     def test_inversion_noisy_sf_huss(self):
         old_model_sf = cfg.PARAMS['use_shape_factor_for_fluxbasedmodel']
         old_inversion_sf = cfg.PARAMS['use_shape_factor_for_inversion']
@@ -1866,6 +1867,7 @@ def test_inversion_tributary_sf_adhikari(self):
         cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = old_model_sf
         cfg.PARAMS['use_shape_factor_for_inversion'] = old_inversion_sf
 
+    @is_slow
     def test_inversion_tributary_sf_huss(self):
         old_model_sf = cfg.PARAMS['use_shape_factor_for_fluxbasedmodel']
         old_inversion_sf = cfg.PARAMS['use_shape_factor_for_inversion']

oggm/tests/test_prepro.py

@@ -207,6 +207,36 @@ def test_glacier_masks(self):
             np.testing.assert_allclose(area*10**-6, gdir.rgi_area_km2,
                                        rtol=1e-1)
 
+    def test_simple_glacier_masks(self):
+
+        # The GIS was double checked externally with IDL.
+        hef_file = get_demo_file('Hintereisferner.shp')
+        entity = gpd.GeoDataFrame.from_file(hef_file).iloc[0]
+
+        gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir)
+        gis.define_glacier_region(gdir, entity=entity)
+        gis.simple_glacier_masks(gdir)
+
+        with netCDF4.Dataset(gdir.get_filepath('gridded_data')) as nc:
+            area = np.sum(nc.variables['glacier_mask'][:] * gdir.grid.dx**2)
+            np.testing.assert_allclose(area*10**-6, gdir.rgi_area_km2,
+                                       rtol=1e-1)
+
+            # Check that HEF doesn't "badly" need a divide
+            mask = nc.variables['glacier_mask'][:]
+            ext = nc.variables['glacier_ext'][:]
+            dem = nc.variables['topo'][:]
+            np.testing.assert_allclose(np.max(dem[mask.astype(bool)]),
+                                       np.max(dem[ext.astype(bool)]),
+                                       atol=10)
+
+        df = utils.glacier_characteristics([gdir], path=False)
+        np.testing.assert_allclose(df['dem_max_elev_on_ext'],
+                                   df['dem_max_elev'],
+                                   atol=10)
+
+        assert np.all(df['dem_max_elev'] > df['dem_max_elev_on_ext'])
+
     def test_intersects(self):
 
         hef_file = get_demo_file('Hintereisferner_RGI5.shp')
@@ -547,6 +577,12 @@ def test_flowlines(self):
         np.testing.assert_allclose(dis * 0 + cfg.PARAMS['flowline_dx'], dis,
                                    rtol=0.01)
 
+        d = gdir.get_diagnostics()
+        assert d['perc_invalid_flowline'] > 0.1
+
+        df = utils.glacier_characteristics([gdir], path=False)
+        assert np.all(df['perc_invalid_flowline'] > 0.1)
+
     def test_geom_width(self):
 
         hef_file = get_demo_file('Hintereisferner_RGI5.shp')