More workflow tools (#264)

oggm/cfg.py

@@ -379,6 +379,13 @@ def oggm_static_paths():
         config['dl_cache_readonly'] = ro
     if os.environ.get('OGGM_DOWNLOAD_CACHE') is not None:
         config['dl_cache_dir'] = os.environ.get('OGGM_DOWNLOAD_CACHE')
+    if os.environ.get('OGGM_EXTRACT_DIR') is not None:
+        # This is for the directories where OGGM needs to extract things
+        # On the cluster it might be useful to do it on a fast disc
+        edir = os.path.abspath(os.environ.get('OGGM_EXTRACT_DIR'))
+        config['tmp_dir'] = os.path.join(edir, 'tmp')
+        config['cru_dir'] = os.path.join(edir, 'cru')
+        config['rgi_dir'] = os.path.join(edir, 'rgi')
 
     if not config['dl_cache_dir']:
         raise RuntimeError('At the very least, the "dl_cache_dir" entry '

oggm/core/preprocessing/gis.py

@@ -99,8 +99,10 @@ def _check_geometry(geometry):
             if parts[0].contains(p):
                 interiors.append(p.exterior)
             else:
-                # This should not happen See that we have a small geom here
-                assert p.area < 1e-9
+                # This should not happen. Check that we have a small geom here
+                if p.area > 1e-4:
+                    log.warning('warning while correcting geometry. Area was: '
+                                '{} but it should be smaller.'.format(p.area))
         geometry = shpg.Polygon(exterior, interiors)
 
     assert 'Polygon' in geometry.type

oggm/sandbox/gmd_paper/plot_hef_dynamics.py

@@ -0,0 +1,82 @@
+import os
+import geopandas as gpd
+import numpy as np
+import oggm
+from oggm import cfg, tasks, graphics
+from oggm.utils import get_demo_file
+import matplotlib.pyplot as plt
+import xarray as xr
+import scipy.optimize as optimization
+from oggm.sandbox.gmd_paper import PLOT_DIR
+from oggm.core.preprocessing.climate import (t_star_from_refmb,
+                                             local_mustar_apparent_mb)
+from oggm.core.preprocessing.inversion import (mass_conservation_inversion)
+from oggm.core.models.flowline import (FluxBasedModel)
+from oggm.core.models.massbalance import (RandomMassBalanceModel)
+
+cfg.initialize()
+
+# Don't use divides for now, or intersects
+cfg.set_divides_db()
+cfg.set_intersects_db()
+
+cfg.PATHS['dem_file'] = get_demo_file('hef_srtm.tif')
+cfg.PARAMS['border'] = 60
+cfg.PARAMS['auto_skip_task'] = True
+
+base_dir = os.path.join(os.path.expanduser('~/tmp'), 'OGGM_GMD', 'dynamics')
+entity = gpd.read_file(get_demo_file('Hintereisferner_RGI5.shp')).iloc[0]
+gdir = oggm.GlacierDirectory(entity, base_dir=base_dir)
+
+tasks.define_glacier_region(gdir, entity=entity)
+tasks.glacier_masks(gdir)
+tasks.compute_centerlines(gdir)
+tasks.compute_downstream_lines(gdir)
+tasks.initialize_flowlines(gdir)
+tasks.compute_downstream_bedshape(gdir)
+tasks.catchment_area(gdir)
+tasks.catchment_intersections(gdir)
+tasks.catchment_width_geom(gdir)
+tasks.catchment_width_correction(gdir)
+tasks.process_cru_data(gdir)
+tasks.mu_candidates(gdir, div_id=0)
+tasks.compute_ref_t_stars([gdir])
+tasks.distribute_t_stars([gdir])
+tasks.prepare_for_inversion(gdir)
+tasks.volume_inversion(gdir, use_cfg_params={'glen_a':cfg.A, 'fs':0})
+tasks.filter_inversion_output(gdir)
+tasks.init_present_time_glacier(gdir)
+
+
+reset = False
+tasks.random_glacier_evolution(gdir, nyears=500, bias=0, seed=0,
+                               filesuffix='seed0', reset=reset)
+
+tasks.random_glacier_evolution(gdir, nyears=500, bias=0, seed=1,
+                               filesuffix='seed1', reset=reset)
+
+tasks.random_glacier_evolution(gdir, nyears=500, bias=0, seed=2,
+                               filesuffix='seed2', reset=reset)
+
+f = gdir.get_filepath('model_diagnostics', filesuffix='seed0')
+ds1 = xr.open_dataset(f)
+f = gdir.get_filepath('model_diagnostics', filesuffix='seed1')
+ds2 = xr.open_dataset(f)
+f = gdir.get_filepath('model_diagnostics', filesuffix='seed2')
+ds3 = xr.open_dataset(f)
+
+# ds1.volume_m3.plot()
+# ds2.volume_m3.plot()
+# ds3.volume_m3.plot()
+# plt.show()
+#
+# ds1.area_m2.plot()
+# ds2.area_m2.plot()
+# ds3.area_m2.plot()
+# plt.show()
+
+ds1.length_m.rolling(time=5, center=True).mean().plot()
+ds2.length_m.rolling(time=5, center=True).mean().plot()
+ds3.length_m.rolling(time=5, center=True).mean().plot()
+plt.show()
+

oggm/sandbox/gmd_paper/plot_inversion_sensi.py

@@ -99,7 +99,5 @@
 ax.set_ylabel('Thickness [m]')
 ax.legend(loc=3)
 
-
-plt.show()
-# plt.tight_layout()
-# plt.savefig(PLOT_DIR + 'inversions_sensi.pdf', dpi=150, bbox_inches='tight')
+plt.tight_layout()
+plt.savefig(PLOT_DIR + 'inversions_sensi.pdf', dpi=150, bbox_inches='tight')

oggm/sandbox/gmd_paper/plot_workflow.py

@@ -65,10 +65,11 @@ def to_optimize(x):
 # run
 tasks.init_present_time_glacier(gdir)
 
-mb = RandomMassBalanceModel(gdir, seed=0)
+mb = RandomMassBalanceModel(gdir, seed=1)
 fls = gdir.read_pickle('model_flowlines')
 model = FluxBasedModel(fls, mb_model=mb, y0=0, glen_a=glen_a)
-model.run_until(100)
+
+model.run_until(150)
 
 f = plt.figure(figsize=(10, 12))
 from mpl_toolkits.axes_grid1 import ImageGrid

oggm/sandbox/rgi_tools/scripts.py

@@ -0,0 +1,18 @@
+"""This script reads all the RGI files and computes intersects out of them."""
+import multiprocessing as mp
+import os
+from glob import glob
+
+from oggm import cfg
+from oggm.utils import get_rgi_dir
+
+# Download RGI files
+cfg.initialize()
+rgi_dir = get_rgi_dir()
+rgi_shps = list(glob(os.path.join(rgi_dir, "*", '*_rgi50_*.shp')))
+rgi_shps = [r for r in rgi_shps if 'Regions' not in r]
+
+from oggm.sandbox.rgi_tools.tools import compute_intersects
+
+with mp.Pool() as p:
+    p.map(compute_intersects, rgi_shps, chunksize=1)

oggm/sandbox/rgi_tools/tools.py

@@ -0,0 +1,87 @@
+"""This script reads all the RGI files and computes intersects out of them."""
+import os
+import geopandas as gpd
+from shapely.ops import linemerge
+import shapely.geometry as shpg
+from oggm.utils import haversine, mkdir
+from oggm.core.preprocessing.gis import _check_geometry
+from salem import wgs84
+import time
+
+OUTDIR = '/home/mowglie/disk/Data/GIS/SHAPES/RGI/RGI_V5_Intersects'
+
+
+def compute_intersects(rgi_shp):
+    """Processes the rgi file and writes the interscects to out_path"""
+
+    out_path = os.path.basename(rgi_shp)
+    odir = os.path.basename(os.path.dirname(rgi_shp))
+    odir = os.path.join(OUTDIR, odir)
+    mkdir(odir)
+    out_path = os.path.join(odir, 'intersects_' + out_path)
+
+    print('Start ' + os.path.basename(rgi_shp) + ' ...')
+    start_time = time.time()
+
+    gdf = gpd.read_file(rgi_shp)
+
+    # clean geometries like OGGM does
+    gdf['geometry'] = [_check_geometry(g) for g in gdf.geometry]
+
+    out_cols = ['RGIId_1', 'RGIId_2', 'geometry']
+    out = gpd.GeoDataFrame(columns=out_cols)
+
+    for i, major in gdf.iterrows():
+
+        # Exterior only
+        major_poly = major.geometry.exterior
+
+        # sort by distance to the current glacier
+        gdf['dis'] = haversine(major.CenLon, major.CenLat,
+                               gdf.CenLon, gdf.CenLat)
+        gdfs = gdf.sort_values(by='dis').iloc[1:]
+
+        # Keep glaciers in which intersect
+        gdfs = gdfs.loc[gdfs.dis < 200000]
+        gdfs = gdfs.loc[gdfs.intersects(major_poly)]
+
+        for i, neighbor in gdfs.iterrows():
+
+            if neighbor.RGIId in out.RGIId_1 or neighbor.RGIId in out.RGIId_2:
+                continue
+
+            # Exterior only
+            # Buffer is needed for numerical reasons
+            neighbor_poly = neighbor.geometry.exterior.buffer(0.0001)
+
+            # Go
+            try:
+                mult_intersect = major_poly.intersection(neighbor_poly)
+            except:
+                continue
+
+            if isinstance(mult_intersect, shpg.Point):
+                continue
+            if isinstance(mult_intersect, shpg.linestring.LineString):
+                mult_intersect = [mult_intersect]
+            if len(mult_intersect) == 0:
+                continue
+            mult_intersect = [m for m in mult_intersect if
+                              not isinstance(m, shpg.Point)]
+            if len(mult_intersect) == 0:
+                continue
+            mult_intersect = linemerge(mult_intersect)
+            if isinstance(mult_intersect, shpg.linestring.LineString):
+                mult_intersect = [mult_intersect]
+            for line in mult_intersect:
+                assert isinstance(line, shpg.linestring.LineString)
+                line = gpd.GeoDataFrame([[major.RGIId, neighbor.RGIId, line]],
+                                        columns=out_cols)
+                out = out.append(line)
+
+    out.crs = wgs84.srs
+    out.to_file(out_path)
+
+    print(os.path.basename(rgi_shp) +
+          ' took {0:.2f} seconds'.format(time.time() - start_time))
+    return

oggm/tests/test_utils.py

@@ -412,6 +412,19 @@ def test_download_rgi(self):
 
         cfg.PATHS['rgi_dir'] = tmp
 
+    @is_download
+    def test_download_rgi_intersects(self):
+
+        tmp = cfg.PATHS['rgi_dir']
+        cfg.PATHS['rgi_dir'] = os.path.join(TEST_DIR, 'rgi_extract')
+
+        of = utils.get_rgi_intersects_dir()
+        of = os.path.join(of, '01_rgi50_Alaska',
+                          'intersects_01_rgi50_Alaska.shp')
+        self.assertTrue(os.path.exists(of))
+
+        cfg.PATHS['rgi_dir'] = tmp
+
     @is_download
     def test_download_dem3_viewpano(self):
 

oggm/utils.py

@@ -1366,10 +1366,10 @@ def _get_rgi_dir_unlocked():
         raise ValueError('The RGI data directory has to be'
                          'specified explicitly.')
     rgi_dir = os.path.abspath(os.path.expanduser(rgi_dir))
+    rgi_dir = os.path.join(rgi_dir, 'RGIV5')
     mkdir(rgi_dir)
 
-    bname = 'rgi50.zip'
-    dfile = 'http://www.glims.org/RGI/rgi50_files/' + bname
+    dfile = 'http://www.glims.org/RGI/rgi50_files/rgi50.zip'
     test_file = os.path.join(rgi_dir, '000_rgi50_manifest.txt')
 
     if not os.path.exists(test_file):
@@ -1392,6 +1392,46 @@ def _get_rgi_dir_unlocked():
     return rgi_dir
 
 
+def get_rgi_intersects_dir():
+    """Returns a path to the RGI directory containing the intersects.
+
+    If the files are not present, download them.
+
+    Returns
+    -------
+    path to the directory
+    """
+
+    with _get_download_lock():
+        return _get_rgi_intersects_dir_unlocked()
+
+
+def _get_rgi_intersects_dir_unlocked():
+
+    rgi_dir = cfg.PATHS['rgi_dir']
+
+    # Be sure the user gave a sensible path to the RGI dir
+    if not rgi_dir:
+        raise ValueError('The RGI data directory has to be'
+                         'specified explicitly.')
+
+    rgi_dir = os.path.abspath(os.path.expanduser(rgi_dir))
+    mkdir(rgi_dir)
+
+    dfile = ('https://dl.dropboxusercontent.com/u/20930277/OGGM_Public/' +
+             'RGI_V5_Intersects.zip')
+    test_file = os.path.join(rgi_dir, 'RGI_V5_Intersects',
+                             'Intersects_OGGM_Manifest.txt')
+
+    if not os.path.exists(test_file):
+        # if not there download it
+        ofile = file_downloader(dfile)
+        # Extract root
+        with zipfile.ZipFile(ofile) as zf:
+            zf.extractall(rgi_dir)
+    return os.path.join(rgi_dir, 'RGI_V5_Intersects')
+
+
 def get_cru_file(var=None):
     """Returns a path to the desired CRU TS file.