Merge 87b0012 into f8a1745

oggm/cli/prepro_levels.py

@@ -360,7 +360,14 @@ def _time_log():
             workflow.execute_entity_task(gis.rasterio_glacier_mask,
                                          gdirs, source='ALL')
 
-            # Compress all in output directory
+            # Glacier stats
+            sum_dir = os.path.join(output_base_dir, 'L1', 'summary')
+            utils.mkdir(sum_dir)
+            opath = os.path.join(sum_dir, 'glacier_statistics_{}.csv'.format(rgi_reg))
+            utils.compile_glacier_statistics(gdirs, path=opath)
+
+            # L1 OK - compress all in output directory
+            log.workflow('L1 done. Writing to tar...')
             level_base_dir = os.path.join(output_base_dir, 'L1')
             workflow.execute_entity_task(utils.gdir_to_tar, gdirs, delete=False,
                                          base_dir=level_base_dir)

oggm/core/gis.py

@@ -952,16 +952,26 @@ def simple_glacier_masks(gdir, write_hypsometry=False):
     slope = np.arctan(np.sqrt(sx ** 2 + sy ** 2))
     avg_slope = np.rad2deg(np.mean(slope[glacier_mask]))
 
+    sec_bins = -22.5 + 45 * np.arange(9)
+    aspect_for_bin = aspect
+    if aspect_for_bin >= sec_bins[-1]:
+        aspect_for_bin -= 360
+    aspect_sec = np.digitize(aspect_for_bin, sec_bins)
+    dx2 =gdir.grid.dx**2 * 1e-6
     # write
     df = pd.DataFrame()
-    df['RGIId'] = [gdir.rgi_id]
-    df['GLIMSId'] = [gdir.glims_id]
-    df['Zmin'] = [dem_on_ice.min()]
-    df['Zmax'] = [dem_on_ice.max()]
-    df['Zmed'] = [np.median(dem_on_ice)]
-    df['Area'] = [gdir.rgi_area_km2]
-    df['Slope'] = [avg_slope]
-    df['Aspect'] = [aspect]
+    df['rgi_id'] = [gdir.rgi_id]
+    df['area_km2'] = [gdir.rgi_area_km2]
+    df['area_grid_km2'] = [glacier_mask.sum() * dx2]
+    df['valid_dem_area_km2'] = [np.isfinite(dem_on_ice).sum() * dx2]
+    df['zmin_m'] = [np.nanmin(dem_on_ice)]
+    df['zmax_m'] = [np.nanmax(dem_on_ice)]
+    df['zmed_m'] = [np.nanmedian(dem_on_ice)]
+    df['zmean_m'] = [np.nanmean(dem_on_ice)]
+    df['slope_deg'] = [avg_slope]
+    df['aspect_deg'] = [aspect]
+    df['aspect_sec'] = [aspect_sec]
+    df['dem_source'] = [gdir.get_diagnostics()['dem_source']]
     for b, bs in zip(hi, (bins[1:] + bins[:-1])/2):
         df['{}'.format(np.round(bs).astype(int))] = [b]
     df.to_csv(gdir.get_filepath('hypsometry'), index=False)
@@ -1454,7 +1464,7 @@ def proj(x, y):
     if not glacier_poly_hr.is_valid:
         raise RuntimeError('This glacier geometry is not valid.')
 
-    # Rounded geometry to nearest nearest pix
+    # Rounded geometry to nearest pix
     # I can not use _polyg
     # glacier_poly_pix = _polygon_to_pix(glacier_poly_hr)
     def project(x, y):

oggm/shop/rgitopo.py

@@ -90,29 +90,27 @@ def select_dem_from_dir(gdir, dem_source=None, keep_dem_folders=False):
     gdir : :py:class:`oggm.GlacierDirectory`
         the glacier directory
     dem_source : str
-        the source to pick from (default: NASADEM and COPDEM90)
+        the source to pick from
     keep_dem_folders : bool
         the default is to delete the other DEM directories to save space.
         Set this to True to prevent that (e.g. for sensitivity tests)
     """
 
     # Start by deleting noise
-    for fn in ['log.txt', 'diagnostics.json']:
-        fp = os.path.join(gdir.dir, fn)
-        if os.path.exists(fp):
-            os.remove(fp)
+    for fn in os.listdir(gdir.dir):
+        if fn in ['glacier_mask.tif', 'glacier_grid.json',
+                  'outlines.tar.gz', 'intersects.tar.gz']:
+            continue
+        fn = os.path.join(gdir.dir, fn)
+        if os.path.isfile(fn):
+            os.remove(fn)
 
     sources = [f.name for f in os.scandir(gdir.dir) if f.is_dir()
                and not f.name.startswith('.')]
-
-    if dem_source is None:
-        for s in ['NASADEM', 'COPDEM90', 'GIMP', 'REMA', 'TANDEM', 'MAPZEN']:
-            if s in sources:
-                dem_source = s
-
     if dem_source not in sources:
         raise RuntimeError('source {} not in folder'.format(dem_source))
 
+    gdir.add_to_diagnostics('dem_source', dem_source)
     shutil.copyfile(os.path.join(gdir.dir, dem_source, 'dem.tif'),
                     gdir.get_filepath('dem'))
     shutil.copyfile(os.path.join(gdir.dir, dem_source, 'dem_source.txt'),
@@ -123,7 +121,11 @@ def select_dem_from_dir(gdir, dem_source=None, keep_dem_folders=False):
             shutil.rmtree(os.path.join(gdir.dir, source))
 
 
-@utils.entity_task(log, writes=[])
+def _fallback_dem_quality_check(gdir):
+    return dict(rgi_id=gdir.rgi_id)
+
+
+@utils.entity_task(log, writes=[], fallback=_fallback_dem_quality_check)
 def dem_quality_check(gdir):
     """Run a simple quality check on the rgitopo DEMs
 
@@ -146,27 +148,18 @@ def dem_quality_check(gdir):
     sources = [f.name for f in os.scandir(gdir.dir) if f.is_dir()
                and not f.name.startswith('.')]
 
-    out = {}
+    out = dict(rgi_id=gdir.rgi_id)
     for s in sources:
         try:
             with rasterio.open(os.path.join(gdir.dir, s, 'dem.tif'), 'r',
                                driver='GTiff') as ds:
                 topo = ds.read(1).astype(rasterio.float32)
                 topo[topo <= -999.] = np.NaN
                 topo[ds.read_masks(1) == 0] = np.NaN
-
             valid_mask = np.isfinite(topo) & mask
             if np.all(~valid_mask):
                 continue
-
-            z_on_glacier = topo[valid_mask]
-            zmax, zmin = np.nanmax(z_on_glacier), np.nanmin(z_on_glacier)
-            if zmin == zmax:
-                continue
-
             out[s] = valid_mask.sum() / area
-
         except BaseException:
             pass
-
     return out

oggm/tests/test_prepro.py

@@ -321,11 +321,11 @@ def test_simple_glacier_masks(self):
 
         dfh = pd.read_csv(gdir.get_filepath('hypsometry'))
 
-        np.testing.assert_allclose(dfh['Slope'], entity.Slope, atol=0.5)
-        np.testing.assert_allclose(dfh['Aspect'], entity.Aspect, atol=5)
-        np.testing.assert_allclose(dfh['Zmed'], entity.Zmed, atol=20)
-        np.testing.assert_allclose(dfh['Zmax'], entity.Zmax, atol=20)
-        np.testing.assert_allclose(dfh['Zmin'], entity.Zmin, atol=20)
+        np.testing.assert_allclose(dfh['slope_deg'], entity.Slope, atol=0.5)
+        np.testing.assert_allclose(dfh['aspect_deg'], entity.Aspect, atol=5)
+        np.testing.assert_allclose(dfh['zmed_m'], entity.Zmed, atol=20)
+        np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=20)
+        np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=20)
 
         bins = []
         for c in dfh.columns:
@@ -369,11 +369,11 @@ def test_glacier_masks_other_glacier(self):
                         os.path.join(self.testdir, 'simple_masks.nc'))
 
         dfh = pd.read_csv(gdir.get_filepath('hypsometry'))
-        np.testing.assert_allclose(dfh['Slope'], entity.Slope, atol=1)
-        np.testing.assert_allclose(dfh['Aspect'], entity.Aspect, atol=10)
-        np.testing.assert_allclose(dfh['Zmed'], entity.Zmed, atol=20)
-        np.testing.assert_allclose(dfh['Zmax'], entity.Zmax, atol=20)
-        np.testing.assert_allclose(dfh['Zmin'], entity.Zmin, atol=20)
+        np.testing.assert_allclose(dfh['slope_deg'], entity.Slope, atol=1)
+        np.testing.assert_allclose(dfh['aspect_deg'], entity.Aspect, atol=10)
+        np.testing.assert_allclose(dfh['zmed_m'], entity.Zmed, atol=20)
+        np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=20)
+        np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=20)
 
     @pytest.mark.skipif((Version(rasterio.__version__) <
                          Version('1.0')),

oggm/tests/test_shop.py

@@ -214,7 +214,8 @@ def test_from_dem(self, class_case_dir, monkeypatch):
                                                  'emen.de/~oggm/test_gdirs/dem'
                                                  's_v1/default/')
 
-        gd = rgitopo.init_glacier_directories_from_rgitopo(['RGI60-09.01004'])
+        gd = rgitopo.init_glacier_directories_from_rgitopo(['RGI60-09.01004'],
+                                                           dem_source='COPDEM')
         gd = gd[0]
 
         assert gd.has_file('dem')
@@ -237,10 +238,11 @@ def test_qc(self, class_case_dir, monkeypatch):
                                                  's_v1/default/')
 
         gd = rgitopo.init_glacier_directories_from_rgitopo(['RGI60-09.01004'],
+                                                           dem_source='COPDEM',
                                                            keep_dem_folders=True)
         out = rgitopo.dem_quality_check(gd[0])
         assert len(out) > 5
-        assert np.sum(list(out.values())) > 5
+        assert pd.Series(out).iloc[1:].sum() > 5
 
 
 class Test_w5e5:

oggm/utils/_downloads.py

@@ -1297,7 +1297,11 @@ def _get_prepro_gdir_unlocked(rgi_version, rgi_id, border, prepro_level,
 
     url = get_prepro_base_url(rgi_version=rgi_version, border=border,
                               prepro_level=prepro_level, base_url=base_url)
-    url += '{}/{}.tar' .format(rgi_id[:8], rgi_id[:11])
+    if len(rgi_id) == 23:
+        # RGI7
+        url += '{}/{}.tar'.format(rgi_id[:17], rgi_id[:20])
+    else:
+        url += '{}/{}.tar'.format(rgi_id[:8], rgi_id[:11])
     tar_base = file_downloader(url)
     if tar_base is None:
         raise RuntimeError('Could not find file at ' + url)

oggm/utils/_funcs.py

@@ -22,6 +22,7 @@
 from scipy.interpolate import interp1d
 import shapely.geometry as shpg
 from shapely.ops import linemerge
+from shapely.validation import make_valid
 
 # Optional libs
 try:
@@ -519,9 +520,29 @@ def polygon_intersections(gdf):
     return out
 
 
+def recursive_valid_polygons(geoms, crs):
+    """Given a list of shapely geometries, makes sure all geometries are valid
+
+    All will be valid polygons of area > 10000m2"""
+    new_geoms = []
+    for geom in geoms:
+        new_geom = make_valid(geom)
+        try:
+            new_geoms.extend(recursive_valid_polygons(list(new_geom.geoms), crs))
+        except AttributeError:
+            new_s = gpd.GeoSeries(new_geom)
+            new_s.crs = crs
+            if new_s.to_crs({'proj': 'cea'}).area.iloc[0] >= 10000:
+                new_geoms.append(new_geom)
+    assert np.all([type(geom) == shpg.Polygon for geom in new_geoms])
+    return new_geoms
+
+
 def multipolygon_to_polygon(geometry, gdir=None):
     """Sometimes an RGI geometry is a multipolygon: this should not happen.
 
+    It was vor vey old versions, pretty sure this is not needed anymore.
+
     Parameters
     ----------
     geometry : shpg.Polygon or shpg.MultiPolygon

oggm/utils/_workflow.py

@@ -60,7 +60,8 @@
 from oggm import __version__
 from oggm.utils._funcs import (calendardate_to_hydrodate, date_to_floatyear,
                                tolist, filter_rgi_name, parse_rgi_meta,
-                               haversine, multipolygon_to_polygon, clip_scalar)
+                               haversine, multipolygon_to_polygon,
+                               recursive_valid_polygons)
 from oggm.utils._downloads import (get_demo_file, get_wgms_files,
                                    get_rgi_glacier_entities)
 from oggm import cfg
@@ -1768,6 +1769,64 @@ def compile_glacier_statistics(gdirs, filesuffix='', path=True,
     return out
 
 
+@entity_task(log)
+def read_glacier_hypsometry(gdir):
+    """Utility function to read the glacier hypsometry in the folder.
+
+    Parameters
+    ----------
+    gdir :  :py:class:`oggm.GlacierDirectory` object
+        the glacier directory to process
+
+    Returns
+    -------
+    the dataframe
+    """
+    try:
+        out = pd.read_csv(gdir.get_filepath('hypsometry')).iloc[0]
+    except:
+        out = pd.Series({'rgi_id': gdir.rgi_id})
+    return out
+
+
+@global_task(log)
+def compile_glacier_hypsometry(gdirs, filesuffix='', path=True,
+                               add_column=None):
+    """Gather as much statistics as possible about a list of glaciers.
+
+    It can be used to do result diagnostics and other stuffs. If the data
+    necessary for a statistic is not available (e.g.: flowlines length) it
+    will simply be ignored.
+
+    Parameters
+    ----------
+    gdirs : list of :py:class:`oggm.GlacierDirectory` objects
+        the glacier directories to process
+    filesuffix : str
+        add suffix to output file
+    path : str, bool
+        Set to "True" in order  to store the info in the working directory
+        Set to a path to store the file to your chosen location
+    add_column : tuple
+        if you feel like adding a key - value pair to the compiled dataframe
+    """
+    from oggm.workflow import execute_entity_task
+
+    out_df = execute_entity_task(read_glacier_hypsometry, gdirs)
+
+    out = pd.DataFrame(out_df).set_index('rgi_id')
+    if add_column is not None:
+        out[add_column[0]] = add_column[1]
+    if path:
+        if path is True:
+            out.to_csv(os.path.join(cfg.PATHS['working_dir'],
+                                    ('glacier_hypsometry' +
+                                     filesuffix + '.csv')))
+        else:
+            out.to_csv(path)
+    return out
+
+
 @global_task(log)
 def compile_fixed_geometry_mass_balance(gdirs, filesuffix='',
                                         path=True, csv=False,
@@ -2593,6 +2652,7 @@ def __init__(self, rgi_entity, base_dir=None, reset=False,
         # Do we want to use the RGI center point or ours?
         if cfg.PARAMS['use_rgi_area']:
             try:
+                # RGIv7
                 self.cenlon = float(rgi_entity.cenlon)
                 self.cenlat = float(rgi_entity.cenlat)
             except AttributeError:
@@ -2782,21 +2842,23 @@ def __repr__(self):
         return '\n'.join(summary) + '\n'
 
     def _reproject_and_write_shapefile(self, entity):
-
         # Make a local glacier map
         if cfg.PARAMS['map_proj'] == 'utm':
-            from pyproj.aoi import AreaOfInterest
-            from pyproj.database import query_utm_crs_info
-            utm_crs_list = query_utm_crs_info(
-                datum_name="WGS 84",
-                area_of_interest=AreaOfInterest(
-                    west_lon_degree=self.cenlon,
-                    south_lat_degree=self.cenlat,
-                    east_lon_degree=self.cenlon,
-                    north_lat_degree=self.cenlat,
-                ),
-            )
-            proj4_str = utm_crs_list[0].code
+            if entity.get('utm_zone', False):
+                proj4_str = {'proj': 'utm', 'zone': entity['utm_zone']}
+            else:
+                from pyproj.aoi import AreaOfInterest
+                from pyproj.database import query_utm_crs_info
+                utm_crs_list = query_utm_crs_info(
+                    datum_name="WGS 84",
+                    area_of_interest=AreaOfInterest(
+                        west_lon_degree=self.cenlon,
+                        south_lat_degree=self.cenlat,
+                        east_lon_degree=self.cenlon,
+                        north_lat_degree=self.cenlat,
+                    ),
+                )
+                proj4_str = utm_crs_list[0].code
         elif cfg.PARAMS['map_proj'] == 'tmerc':
             params = dict(name='tmerc', lat_0=0., lon_0=self.cenlon,
                           k=0.9996, x_0=0, y_0=0, datum='WGS84')
@@ -2812,6 +2874,16 @@ def _reproject_and_write_shapefile(self, entity):
         # transform geometry to map
         project = partial(transform_proj, proj_in, proj_out)
         geometry = shp_trafo(project, entity['geometry'])
+        if len(self.rgi_id) == 23 and not geometry.is_valid:
+            # In RGI7 we know that the geometries are valid in entry
+            # So we have to validate them after proj as well
+            # Try buffer first
+            geometry = geometry.buffer(0)
+            if not geometry.is_valid:
+                correct = recursive_valid_polygons([geometry], crs=proj4_str)
+                if len(correct) != 1:
+                    raise RuntimeError('Cant correct this geometry')
+                geometry = correct[0]
         geometry = multipolygon_to_polygon(geometry, gdir=self)
 
         # Save transformed geometry to disk
@@ -2886,7 +2958,7 @@ def rgi_area_km2(self):
         except KeyError:
             # RGI V7
             _area = self.read_shapefile('outlines')['area_km2']
-        return np.round(float(_area.iloc[0]), decimals=3)
+        return float(_area.iloc[0])
 
     @lazy_property
     def intersects_ids(self):