More options for prepro directories (#1664)

oggm/cfg.py

@@ -300,7 +300,7 @@ def _log_param_change(self, key, value):
 BASENAMES['elevation_band_flowline'] = ('elevation_band_flowline.csv', _doc)
 
 
-def set_logging_config(logging_level='INFO', future=False):
+def set_logging_config(logging_level='INFO'):
     """Set the global logger parameters.
 
     Logging levels:
@@ -329,8 +329,6 @@ def set_logging_config(logging_level='INFO', future=False):
         the logging level. See description above for a list of options. Setting
         to `None` is equivalent to `'CRITICAL'`, i.e. no log output will be
         generated.
-    future : bool
-        use the new behavior of logging='WORKFLOW'.
     """
 
     # Add a custom level - just for us
@@ -363,31 +361,12 @@ def workflow(self, message, *args, **kws):
 
     logging_level = logging_level.upper()
 
-    # Deprecation warning
-    if logging_level == 'WORKFLOW' and not future:
-
-        msg = ('In future versions of OGGM, the logging config WORKFLOW '
-               'will no longer print ERROR or WARNING messages, but only high '
-               'level information (i.e. hiding potential errors in your code '
-               'but also avoiding cluttered log files for runs with '
-               'many expected errors, e.g. global runs). If you want to obtain '
-               'a similar logger behavior as before, set '
-               "`logging_level='WARNING'`, which will print high level info "
-               "as well as errors and warnings during the run. If you "
-               "want to use the new behavior and suppress this warning, "
-               "set `logging_level='WORKFLOW'` and `future=True`.")
-        warnings.warn(msg, category=FutureWarning)
-
-        # Set old behavior
-        logging_level = 'WARNING'
-
     logging.basicConfig(format='%(asctime)s: %(name)s: %(message)s',
                         datefmt='%Y-%m-%d %H:%M:%S',
                         level=getattr(logging, logging_level))
 
 
-def initialize_minimal(file=None, logging_level='INFO', params=None,
-                       future=False):
+def initialize_minimal(file=None, logging_level='INFO', params=None):
     """Same as initialise() but without requiring any download of data.
 
     This is useful for "flowline only" OGGM applications
@@ -400,14 +379,12 @@ def initialize_minimal(file=None, logging_level='INFO', params=None,
         set a logging level. See :func:`set_logging_config` for options.
     params : dict
         overrides for specific parameters from the config file
-    future : bool
-        use the new behavior of logging='WORKFLOW'.
     """
     global IS_INITIALIZED
     global PARAMS
     global PATHS
 
-    set_logging_config(logging_level=logging_level, future=future)
+    set_logging_config(logging_level=logging_level)
 
     is_default = False
     if file is None:
@@ -547,6 +524,10 @@ def initialize_minimal(file=None, logging_level='INFO', params=None,
     PARAMS[k] = [str(vk) for vk in cp.as_list(k)]
     k = 'store_fl_diagnostic_variables'
     PARAMS[k] = [str(vk) for vk in cp.as_list(k)]
+    k = 'by_bin_dx'
+    PARAMS[k] = [float(vk) for vk in cp.as_list(k)]
+    k = 'by_bin_bins'
+    PARAMS[k] = [float(vk) for vk in cp.as_list(k)]
 
     # Flowline model
     k = 'glacier_length_method'
@@ -571,12 +552,12 @@ def initialize_minimal(file=None, logging_level='INFO', params=None,
 
     # Delete non-floats
     ltr = ['working_dir', 'dem_file', 'climate_file', 'use_tar_shapefiles',
-           'grid_dx_method', 'compress_climate_netcdf',
+           'grid_dx_method', 'compress_climate_netcdf', 'by_bin_dx',
            'mp_processes', 'use_multiprocessing', 'clip_dem_to_zero',
            'topo_interp', 'use_compression', 'bed_shape', 'continue_on_error',
            'use_multiple_flowlines', 'border', 'use_temp_bias_from_file',
            'mpi_recv_buf_size', 'map_proj', 'evolution_model',
-           'hydro_month_sh', 'hydro_month_nh',
+           'hydro_month_sh', 'hydro_month_nh', 'by_bin_bins',
            'use_intersects', 'filter_min_slope', 'clip_tidewater_border',
            'auto_skip_task', 'ref_mb_valid_window',
            'rgi_version', 'dl_verify', 'use_mp_spawn', 'calving_use_limiter',
@@ -602,7 +583,7 @@ def initialize_minimal(file=None, logging_level='INFO', params=None,
     IS_INITIALIZED = True
 
 
-def initialize(file=None, logging_level='INFO', params=None, future=False):
+def initialize(file=None, logging_level='INFO', params=None):
     """Read the configuration file containing the run's parameters.
 
     This should be the first call, before using any of the other OGGM modules
@@ -616,20 +597,17 @@ def initialize(file=None, logging_level='INFO', params=None, future=False):
         set a logging level. See :func:`set_logging_config` for options.
     params : dict
         overrides for specific parameters from the config file
-    future : bool
-        use the new behavior of logging='WORKFLOW'.
     """
     global PARAMS
     global DATA
 
-    initialize_minimal(file=file, logging_level=logging_level, params=params,
-                       future=future)
+    initialize_minimal(file=file, logging_level=logging_level, params=params)
 
     # Do not spam
     PARAMS.do_log = False
 
     # Make sure we have a proper cache dir
-    from oggm.utils import download_oggm_files, get_demo_file
+    from oggm.utils import download_oggm_files
     download_oggm_files()
 
     # Read in the demo glaciers

oggm/cli/benchmark.py

@@ -70,7 +70,7 @@ def run_benchmark(rgi_version=None, rgi_reg=None, border=None,
     override_params['working_dir'] = working_dir
 
     # Initialize OGGM and set up the run parameters
-    cfg.initialize(logging_level=logging_level, params=override_params, future=True)
+    cfg.initialize(logging_level=logging_level, params=override_params)
 
     # Use multiprocessing?
     cfg.PARAMS['use_multiprocessing'] = True

oggm/cli/prepro_levels.py

@@ -81,9 +81,10 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
                       elev_bands=False, centerlines=False,
                       override_params=None,
                       mb_calibration_strategy='informed_threestep',
+                      select_source_from_dir=None, keep_dem_folders=False,
                       add_consensus_thickness=False, add_itslive_velocity=False,
                       add_millan_thickness=False, add_millan_velocity=False,
-                      add_hugonnet_dhdt=False,
+                      add_hugonnet_dhdt=False, add_glathida=False,
                       start_level=None, start_base_url=None, max_level=5,
                       logging_level='WORKFLOW',
                       dynamic_spinup=False, err_dmdtda_scaling_factor=0.2,
@@ -102,7 +103,8 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
     output_folder : str
         path to the output folder (where to put the preprocessed tar files)
     dem_source : str
-        which DEM source to use: default, SOURCE_NAME or ALL
+        which DEM source to use: default, SOURCE_NAME, STANDARD or ALL
+        "standard" is COPDEM + NASADEM
     working_dir : str
         path to the OGGM working directory
     params_file : str
@@ -132,6 +134,14 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
         - 'informed_threestep' (default)
         - 'melt_temp'
         - 'temp_melt'
+    select_source_from_dir : str
+        if starting from a level 1 "ALL" or "STANDARD" DEM sources directory,
+        select the chosen DEM source here. If you set it to "BY_RES" here,
+        COPDEM will be used and its resolution chosen based on the gdir's
+        map resolution (COPDEM30 for dx < 60 m, COPDEM90 elsewhere).
+    keep_dem_folders : bool
+        if `select_source_from_dir` is used, wether to keep the original
+        DEM folders in or not.
     add_consensus_thickness : bool
         adds (reprojects) the consensus estimates thickness to the glacier
         directories. With elev_bands=True, the data will also be binned.
@@ -147,6 +157,9 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
     add_hugonnet_dhdt : bool
         adds (reprojects) the hugonnet dhdt maps to the glacier
         directories. With elev_bands=True, the data will also be binned.
+    add_glathida : bool
+        adds (reprojects) the glathida thickness data to the glacier
+        directories. Data points are stored as csv files.
     start_level : int
         the pre-processed level to start from (default is to start from
         scratch). If set, you'll need to indicate start_base_url as well.
@@ -238,8 +251,7 @@ def _time_log():
 
     # Initialize OGGM and set up the run parameters
     cfg.initialize(file=params_file, params=override_params,
-                   logging_level=logging_level,
-                   future=True)
+                   logging_level=logging_level)
 
     # Prepare the download of climate file to be shared across processes
     # TODO
@@ -348,10 +360,16 @@ def _time_log():
             cfg.PATHS['dem_file'] = test_topofile
 
         # Which DEM source?
-        if dem_source.upper() == 'ALL':
+        if dem_source.upper() in ['ALL', 'STANDARD']:
             # This is the complex one, just do the job and leave
-            log.workflow('Running prepro on ALL sources')
-            for i, s in enumerate(utils.DEM_SOURCES):
+
+            if dem_source.upper() == 'ALL':
+                sources = utils.DEM_SOURCES
+            if dem_source.upper() == 'STANDARD':
+                sources = ['COPDEM30', 'COPDEM90', 'NASADEM']
+
+            log.workflow('Running prepro on several sources')
+            for i, s in enumerate(sources):
                 rs = i == 0
                 log.workflow('Running prepro on sources: {}'.format(s))
                 gdirs = workflow.init_glacier_directories(rgidf, reset=rs,
@@ -421,6 +439,16 @@ def _time_log():
                      'N elev bands type: {}.'
                      ''.format(len(gdirs_cent), len(gdirs_band)))
 
+        # If we are coming from a multi-dem setup, let's select it from there
+        if select_source_from_dir is not None:
+            from oggm.shop.rgitopo import select_dem_from_dir
+            workflow.execute_entity_task(select_dem_from_dir, gdirs_band,
+                                         dem_source=select_source_from_dir,
+                                         keep_dem_folders=keep_dem_folders)
+            workflow.execute_entity_task(select_dem_from_dir, gdirs_cent,
+                                         dem_source=select_source_from_dir,
+                                         keep_dem_folders=keep_dem_folders)
+
         # HH2015 method
         workflow.execute_entity_task(tasks.simple_glacier_masks, gdirs_band)
 
@@ -448,6 +476,9 @@ def _time_log():
             from oggm.shop.hugonnet_maps import hugonnet_to_gdir
             workflow.execute_entity_task(hugonnet_to_gdir, gdirs)
             bin_variables.append('hugonnet_dhdt')
+        if add_glathida:
+            from oggm.shop.glathida import glathida_to_gdir
+            workflow.execute_entity_task(glathida_to_gdir, gdirs)
 
         if bin_variables and gdirs_band:
             workflow.execute_entity_task(tasks.elevation_band_flowline,
@@ -503,14 +534,18 @@ def _time_log():
             from oggm.shop.millan22 import compile_millan_statistics
             opath = os.path.join(sum_dir, 'millan_statistics_{}.csv'.format(rgi_reg))
             compile_millan_statistics(gdirs, path=opath)
-        if add_hugonnet_dhdt:
-            from oggm.shop.hugonnet_maps import compile_hugonnet_statistics
-            opath = os.path.join(sum_dir, 'hugonnet_statistics_{}.csv'.format(rgi_reg))
-            compile_hugonnet_statistics(gdirs, path=opath)
         if add_consensus_thickness:
             from oggm.shop.bedtopo import compile_consensus_statistics
             opath = os.path.join(sum_dir, 'consensus_statistics_{}.csv'.format(rgi_reg))
             compile_consensus_statistics(gdirs, path=opath)
+        if add_hugonnet_dhdt:
+            from oggm.shop.hugonnet_maps import compile_hugonnet_statistics
+            opath = os.path.join(sum_dir, 'hugonnet_statistics_{}.csv'.format(rgi_reg))
+            compile_hugonnet_statistics(gdirs, path=opath)
+        if add_glathida:
+            from oggm.shop.glathida import compile_glathida_statistics
+            opath = os.path.join(sum_dir, 'glathida_statistics_{}.csv'.format(rgi_reg))
+            compile_glathida_statistics(gdirs, path=opath)
 
         # And for level 2: shapes
         if len(gdirs_cent) > 0:
@@ -808,6 +843,16 @@ def parse_args(args):
                              'compatible with level 1 folders, after which '
                              'the processing will stop. The default is to use '
                              'the default OGGM DEM.')
+    parser.add_argument('--select-source-from-dir', type=str,
+                        default=None,
+                        help='if starting from a level 1 "ALL" or "STANDARD" DEM '
+                        'sources directory, select the chosen DEM source here. '
+                        'If you set it to "BY_RES" here, COPDEM will be used and '
+                        'its resolution chosen based on the gdirs map resolution '
+                        '(COPDEM30 for dx < 60 m, COPDEM90 elsewhere).')
+    parser.add_argument('--keep-dem-folders', nargs='?', const=True, default=False,
+                        help='if `select_source_from_dir` is used, wether to keep '
+                        'the original DEM folders in or not.')
     parser.add_argument('--add-consensus-thickness', nargs='?', const=True, default=False,
                         help='adds (reprojects) the consensus thickness '
                              'estimates to the glacier directories. '
@@ -833,6 +878,10 @@ def parse_args(args):
                              'maps to the glacier directories. '
                              'With --elev-bands, the data will also be '
                              'binned.')
+    parser.add_argument('--add-glathida', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the glathida point thickness '
+                             'observations to the glacier directories. '
+                             'The data points are stored as csv.')
     parser.add_argument('--demo', nargs='?', const=True, default=False,
                         help='if you want to run the prepro for the '
                              'list of demo glaciers.')
@@ -913,11 +962,14 @@ def parse_args(args):
                 logging_level=args.logging_level,
                 elev_bands=args.elev_bands,
                 centerlines=args.centerlines,
+                select_source_from_dir=args.select_source_from_dir,
+                keep_dem_folders=args.keep_dem_folders,
                 add_consensus_thickness=args.add_consensus_thickness,
                 add_millan_thickness=args.add_millan_thickness,
                 add_itslive_velocity=args.add_itslive_velocity,
                 add_millan_velocity=args.add_millan_velocity,
                 add_hugonnet_dhdt=args.add_hugonnet_dhdt,
+                add_glathida=args.add_glathida,
                 dynamic_spinup=dynamic_spinup,
                 err_dmdtda_scaling_factor=args.err_dmdtda_scaling_factor,
                 dynamic_spinup_start_year=args.dynamic_spinup_start_year,

oggm/core/centerlines.py

@@ -2326,12 +2326,15 @@ def elevation_band_flowline(gdir, bin_variables=None, preserve_totals=True):
         df = df.dropna(how='all', subset=bin_variables)
 
     # Check for binned vars
-    for var, data, in_total, do_p in zip(bin_variables, out_vars, out_totals,
-                                         preserve_totals):
-        if do_p:
-            out_total = np.nansum(df[var] * df['area'])
-            if out_total > 0:
-                df[var] *= in_total / out_total
+    with warnings.catch_warnings():
+        # This can trigger an invalid value
+        warnings.filterwarnings("ignore", category=RuntimeWarning)
+        for var, data, in_total, do_p in zip(bin_variables, out_vars, out_totals,
+                                             preserve_totals):
+            if do_p:
+                out_total = np.nansum(df[var] * df['area'])
+                if out_total > 0:
+                    df[var] *= in_total / out_total
 
     # In OGGM we go from top to bottom
     df = df[::-1]

oggm/core/gis.py

@@ -261,6 +261,13 @@ def glacier_grid_params(gdir):
         dx = np.rint(cfg.PARAMS['d1'] * np.sqrt(area) + cfg.PARAMS['d2'])
     elif dxmethod == 'fixed':
         dx = np.rint(cfg.PARAMS['fixed_dx'])
+    elif dxmethod == 'by_bin':
+        bins = cfg.PARAMS['by_bin_bins']
+        bin_dx = cfg.PARAMS['by_bin_dx']
+        for i, (b1, b2) in enumerate(zip(bins[:-1], bins[1:])):
+            if b1 < area <= b2:
+                dx = np.rint(bin_dx[i])
+                break
     else:
         raise InvalidParamsError('grid_dx_method not supported: {}'
                                  .format(dxmethod))

oggm/params.cfg

@@ -72,6 +72,7 @@ map_proj = 'tmerc'
 
 # Decision on grid spatial resolution for each glacier
 # 'fixed': dx (meters) = fixed_dx
+# 'by_bin': dx (meters) = dx, chosen by bins
 # 'linear':  dx (meters) = d1 * AREA (km) + d2 ; clipped to dmax (e.g.: 5, 10, 200)
 # 'square':  dx (meters) = d1 * sqrt(AREA) (km) + d2 ; clipped to dmax (e.g.: 20, 10, 200)
 
@@ -84,6 +85,12 @@ dmax = 200.
 # Ignored if grid_dx_method != 'fixed'
 fixed_dx = 50.
 
+# Ignored if grid_dx_method != 'by_bin'
+# size_bins should be of len(by_bin_dx)+1 and given in km2.
+# The first and last bin range should be chosen so that all glaciers are included
+by_bin_dx = 25, 50, 100, 200
+by_bin_bins = 0, 8, 80, 300, 1e12
+
 # Which algorithm to use for interpolating the topography to the local grid
 # 'bilinear' or 'cubic'
 topo_interp = cubic