diff --git a/benchmarks/track_model_results.py b/benchmarks/track_model_results.py index febd6fae3..68717f2e8 100644 --- a/benchmarks/track_model_results.py +++ b/benchmarks/track_model_results.py @@ -36,7 +36,7 @@ def setup_cache(self): gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - tasks.define_glacier_region(gdir, entity=entity) + tasks.define_glacier_region(gdir) tasks.glacier_masks(gdir) tasks.compute_centerlines(gdir) tasks.initialize_flowlines(gdir) @@ -172,10 +172,11 @@ def setup_cache(self): rgidf = rgidf.sort_values('Area', ascending=False) # Go - initialize glacier directories - gdirs = workflow.init_glacier_regions(rgidf) + gdirs = workflow.init_glacier_directories(rgidf) # Preprocessing tasks task_list = [ + tasks.define_glacier_region, tasks.glacier_masks, tasks.compute_centerlines, tasks.initialize_flowlines, @@ -276,7 +277,7 @@ def setup_cache(self): entity = gpd.read_file(get_demo_file('01_rgi60_Columbia.shp')).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - tasks.define_glacier_region(gdir, entity=entity) + tasks.define_glacier_region(gdir) tasks.glacier_masks(gdir) tasks.compute_centerlines(gdir) tasks.initialize_flowlines(gdir) diff --git a/docs/_code/prepare_climate.py b/docs/_code/prepare_climate.py index 1624d79f0..99d422f52 100644 --- a/docs/_code/prepare_climate.py +++ b/docs/_code/prepare_climate.py @@ -22,7 +22,7 @@ entity = gpd.read_file(get_demo_file('HEF_MajDivide.shp')).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=base_dir, reset=True) -tasks.define_glacier_region(gdir, entity=entity) +tasks.define_glacier_region(gdir) tasks.glacier_masks(gdir) tasks.compute_centerlines(gdir) tasks.initialize_flowlines(gdir) diff --git a/docs/_code/prepare_flowlines.py b/docs/_code/prepare_flowlines.py index dc282abc9..58ec13921 100644 --- a/docs/_code/prepare_flowlines.py +++ b/docs/_code/prepare_flowlines.py @@ -12,7 +12,7 @@ entity = gpd.read_file(get_demo_file('HEF_MajDivide.shp')).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=base_dir, reset=True) -tasks.define_glacier_region(gdir, entity=entity) +tasks.define_glacier_region(gdir) tasks.glacier_masks(gdir) tasks.compute_centerlines(gdir) tasks.initialize_flowlines(gdir) diff --git a/docs/api.rst b/docs/api.rst index d484279ce..257a58bda 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -25,6 +25,7 @@ Tools to set-up and run OGGM. cfg.set_logging_config cfg.set_intersects_db cfg.reset_working_dir + workflow.init_glacier_directories workflow.init_glacier_regions workflow.execute_entity_task workflow.gis_prepro_tasks @@ -174,9 +175,10 @@ pre-processed state available on the OGGM servers: cfg.initialize() # always initialize before an OGGM task # The working directory is where OGGM will store the run's data cfg.PATHS['working_dir'] = os.path.join(gettempdir(), 'Docs_GlacierDir') - gdirs = workflow.init_glacier_regions('RGI60-11.00897', from_prepro_level=1, - prepro_border=10) - gdir = gdirs[0] # init_glacier_regions always returns a list + gdirs = workflow.init_glacier_directories('RGI60-11.00897', + from_prepro_level=1, + prepro_border=10) + gdir = gdirs[0] # init_glacier_directories always returns a list We just downloaded the minimal input for a glacier directory. The @@ -353,9 +355,10 @@ including the model flowlines. This is achieved by choosing preprocessing level cfg.initialize() # always initialize before an OGGM task # The working directory is where OGGM will store the run's data cfg.PATHS['working_dir'] = os.path.join(gettempdir(), 'Docs_GlacierDir2') - gdirs = workflow.init_glacier_regions('RGI60-11.00897', from_prepro_level=4, - prepro_border=10) - gdir = gdirs[0] # init_glacier_regions always returns a list + gdirs = workflow.init_glacier_directories('RGI60-11.00897', + from_prepro_level=4, + prepro_border=10) + gdir = gdirs[0] # init_glacier_directories always returns a list fls = gdir.read_pickle('model_flowlines') fls diff --git a/docs/input-data.rst b/docs/input-data.rst index 866e068e6..ebcdcca3f 100644 --- a/docs/input-data.rst +++ b/docs/input-data.rst @@ -127,7 +127,7 @@ of the current stable OGGM version. If you want to change these parameters, you'll have to do a full run from scratch using the :ref:`rawdata`. To start from a pre-processed state, simply use the -:py:func:`workflow.init_glacier_regions` function with the +:py:func:`workflow.init_glacier_directories` function with the ``from_prepro_level`` and ``prepro_border`` keyword arguments set to the values of your choice. @@ -178,9 +178,9 @@ Here is an example with the Hintereisferner in the Alps: f, axs = plt.subplots(2, 2, figsize=(8, 6)) for ax, border in zip(np.array(axs).flatten(), [10, 80, 160, 250]): - gdir = workflow.init_glacier_regions('RGI60-11.00897', - from_prepro_level=1, - prepro_border=border) + gdir = workflow.init_glacier_directories('RGI60-11.00897', + from_prepro_level=1, + prepro_border=border) graphics.plot_domain(gdir, ax=ax, title='Border: {}'.format(border), add_colorbar=False, lonlat_contours_kwargs={'add_tick_labels':False}) @@ -290,7 +290,7 @@ EndDate not included For Greenland and Antarctica, OGGM does not take into account the connectivity level between the Glaciers and the Ice sheets. We recommend to the users to think about this before they -run the task: ``workflow.init_glacier_regions()``. +run the task: ``workflow.init_glacier_directories``. .. _Randolph Glacier Inventory (RGI): https://www.glims.org/RGI/ @@ -490,8 +490,9 @@ have access to the timeseries through the glacier directory: .. ipython:: python - gdir = workflow.init_glacier_regions('RGI60-11.00897', from_prepro_level=4, - prepro_border=10)[0] + gdir = workflow.init_glacier_directories('RGI60-11.00897', + from_prepro_level=4, + prepro_border=10)[0] mb = gdir.get_ref_mb_data() @savefig plot_ref_mbdata.png width=100% mb[['ANNUAL_BALANCE']].plot(title='WGMS data: Hintereisferner') diff --git a/docs/practicalities.rst b/docs/practicalities.rst index b6bf3c1b3..7489ed0a3 100644 --- a/docs/practicalities.rst +++ b/docs/practicalities.rst @@ -250,7 +250,7 @@ you have to deal with: :py:func:`utils.gdir_to_tar` and :py:func:`utils.base_dir_to_tar` functions to create compressed, aggregated files of your directories. You can later initialize new directories from these tar files with the `from_tar` - keyword argument in :py:func:`workflow.init_glacier_regions`. + keyword argument in :py:func:`workflow.init_glacier_directories`. Run per RGI region, not globally diff --git a/docs/run_examples/_code/mb_crossval.py b/docs/run_examples/_code/mb_crossval.py index a12b2bda5..357b423da 100644 --- a/docs/run_examples/_code/mb_crossval.py +++ b/docs/run_examples/_code/mb_crossval.py @@ -36,7 +36,7 @@ rgidf = gpd.read_file(os.path.join(WORKING_DIR, 'mb_ref_glaciers.shp')) # Go - initialize glacier directories -gdirs = workflow.init_glacier_regions(rgidf) +gdirs = workflow.init_glacier_directories(rgidf) # Cross-validation file = os.path.join(cfg.PATHS['working_dir'], 'ref_tstars.csv') diff --git a/docs/run_examples/_code/run_errors.py b/docs/run_examples/_code/run_errors.py index de9196283..5de2ae386 100644 --- a/docs/run_examples/_code/run_errors.py +++ b/docs/run_examples/_code/run_errors.py @@ -36,7 +36,7 @@ log.workflow('Number of glaciers: {}'.format(len(rgi_ids))) # Go - get the pre-processed glacier directories -gdirs = workflow.init_glacier_regions(rgi_ids, from_prepro_level=4) +gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=4) # We can step directly to the experiment! # Random climate representative for the recent climate (1985-2015) diff --git a/docs/run_examples/_code/run_inversion.py b/docs/run_examples/_code/run_inversion.py index 976ab4a0d..378f56d88 100644 --- a/docs/run_examples/_code/run_inversion.py +++ b/docs/run_examples/_code/run_inversion.py @@ -40,8 +40,9 @@ # Go - get the pre-processed glacier directories # We start at level 3, because we need all data for the inversion -gdirs = workflow.init_glacier_regions(rgidf, from_prepro_level=3, - prepro_border=10) +gdirs = workflow.init_glacier_directories(rgidf, + from_prepro_level=3, + prepro_border=10) # Default parameters # Deformation: from Cuffey and Patterson 2010 diff --git a/docs/run_examples/_code/run_reference_mb_glaciers.py b/docs/run_examples/_code/run_reference_mb_glaciers.py index fcfc5be99..532256779 100644 --- a/docs/run_examples/_code/run_reference_mb_glaciers.py +++ b/docs/run_examples/_code/run_reference_mb_glaciers.py @@ -70,7 +70,7 @@ # We have to check which of them actually have enough mb data. # Let OGGM do it: -gdirs = workflow.init_glacier_regions(rgidf) +gdirs = workflow.init_glacier_directories(rgidf) # We need to know which period we have data for log.info('Process the climate data...') @@ -96,10 +96,11 @@ rgidf = rgidf.sort_values('Area', ascending=False) # Go - initialize glacier directories -gdirs = workflow.init_glacier_regions(rgidf) +gdirs = workflow.init_glacier_directories(rgidf) # Prepro tasks task_list = [ + tasks.define_glacier_region, tasks.glacier_masks, tasks.compute_centerlines, tasks.initialize_flowlines, diff --git a/docs/run_examples/_code/run_rgi_region.py b/docs/run_examples/_code/run_rgi_region.py index 020fe07f3..2220251f8 100644 --- a/docs/run_examples/_code/run_rgi_region.py +++ b/docs/run_examples/_code/run_rgi_region.py @@ -52,7 +52,7 @@ log.workflow('Number of glaciers: {}'.format(len(rgidf))) # Go - get the pre-processed glacier directories -gdirs = workflow.init_glacier_regions(rgidf, from_prepro_level=4) +gdirs = workflow.init_glacier_directories(rgidf, from_prepro_level=4) # We can step directly to a new experiment! # Random climate representative for the recent climate (1985-2015) diff --git a/docs/run_examples/_code/run_with_spinup.py b/docs/run_examples/_code/run_with_spinup.py index 9ae38b2f2..df5a43f88 100644 --- a/docs/run_examples/_code/run_with_spinup.py +++ b/docs/run_examples/_code/run_with_spinup.py @@ -33,7 +33,8 @@ cfg.PARAMS['border'] = 80 # Go - initialize glacier directories -gdirs = workflow.init_glacier_regions(['RGI60-11.00897'], from_prepro_level=4) +gdirs = workflow.init_glacier_directories(['RGI60-11.00897'], + from_prepro_level=4) # Additional climate file (CESM) cfg.PATHS['cesm_temp_file'] = get_demo_file('cesm.TREFHT.160001-200512' diff --git a/docs/run_examples/run_errors.rst b/docs/run_examples/run_errors.rst index 0a6c90c4d..46af20bee 100644 --- a/docs/run_examples/run_errors.rst +++ b/docs/run_examples/run_errors.rst @@ -25,17 +25,20 @@ applications. If everything went "well", you should see an output similar to:: - 2019-02-21 15:51:39: oggm.cfg: Using configuration file: /home/mowglie/Documents/git/oggm-fork/oggm/params.cfg - 2019-02-21 15:51:39: __main__: Starting OGGM run - 2019-02-21 15:51:39: __main__: Number of glaciers: 3 - 2019-02-21 15:51:39: oggm.workflow: init_glacier_regions from prepro level 4 on 3 glaciers. - 2019-02-21 15:51:39: oggm.workflow: Execute entity task gdir_from_prepro on 3 glaciers - 2019-02-21 15:51:39: oggm.workflow: Multiprocessing: using all available processors (N=8) - 2019-02-21 15:51:40: oggm.workflow: Execute entity task run_random_climate on 3 glaciers - 2019-02-21 15:51:40: oggm.core.flowline: RuntimeError occurred during task run_random_climate on RGI60-11.03295: Need a valid `model_flowlines` file. If you explicitly want to use `inversion_flowlines`, set use_inversion_flowlines=True. - 2019-02-21 15:52:11: oggm.core.flowline: RuntimeError occurred during task run_random_climate on RGI60-11.00897: Glacier exceeds domain boundaries. - 2019-02-21 15:52:11: oggm.workflow: Execute entity task glacier_statistics on 3 glaciers - 2019-02-21 15:52:11: __main__: OGGM is done! Time needed: 0:00:32 + 2020-03-24 18:54:32: oggm.cfg: Using configuration file: /home/mowglie/disk/Dropbox/HomeDocs/git/oggm-fork/oggm/params.cfg + 2020-03-24 18:54:32: oggm.cfg: Multiprocessing switched ON according to the parameter file. + 2020-03-24 18:54:32: oggm.cfg: Multiprocessing: using all available processors (N=8) + 2020-03-24 18:54:32: __main__: Starting OGGM run + 2020-03-24 18:54:32: __main__: Number of glaciers: 3 + 2020-03-24 18:54:32: oggm.workflow: init_glacier_directories from prepro level 4 on 3 glaciers. + 2020-03-24 18:54:32: oggm.workflow: Execute entity task gdir_from_prepro on 3 glaciers + 2020-03-24 18:54:32: oggm.workflow: Initializing multiprocessing pool with N=8 processes. + 2020-03-24 18:55:00: oggm.workflow: Execute entity task run_random_climate on 3 glaciers + 2020-03-24 18:55:00: oggm.core.flowline: RuntimeError occurred during task run_random_climate on RGI60-11.03295: Need a valid `model_flowlines` file. If you explicitly want to use `inversion_flowlines`, set use_inversion_flowlines=True. + 2020-03-24 18:55:10: oggm.core.flowline: RuntimeError occurred during task run_random_climate on RGI60-11.00897: Glacier exceeds domain boundaries, at year: 135.83333333333334 + 2020-03-24 18:55:31: oggm.workflow: Execute entity task glacier_statistics on 3 glaciers + 2020-03-24 18:55:31: __main__: OGGM is done! Time needed: 0:00:59 + Dealing with errors ------------------- diff --git a/docs/run_examples/run_inversion.rst b/docs/run_examples/run_inversion.rst index 594a26f7b..9ed75bf93 100644 --- a/docs/run_examples/run_inversion.rst +++ b/docs/run_examples/run_inversion.rst @@ -31,18 +31,18 @@ Here, we run the inversion algorithm for all glaciers in the Pyrenees. If everything went well, you should see an output similar to:: - 2019-02-16 22:16:19: oggm.cfg: Using configuration file: /home/mowglie/Documents/git/oggm-fork/oggm/params.cfg - 2019-02-16 22:16:20: __main__: Starting OGGM inversion run - 2019-02-16 22:16:20: __main__: Number of glaciers: 35 - 2019-02-16 22:16:20: oggm.workflow: init_glacier_regions from prepro level 3 on 35 glaciers. - 2019-02-16 22:16:20: oggm.workflow: Execute entity task gdir_from_prepro on 35 glaciers - 2019-02-16 22:16:20: oggm.workflow: Multiprocessing: using all available processors (N=8) - 2019-02-16 22:16:21: oggm.workflow: Execute entity task mass_conservation_inversion on 35 glaciers - 2019-02-16 22:16:21: oggm.workflow: Execute entity task filter_inversion_output on 35 glaciers + 2020-03-25 12:31:52: oggm.cfg: Using configuration file: /home/mowglie/disk/Dropbox/HomeDocs/git/oggm-fork/oggm/params.cfg + 2020-03-25 12:31:52: oggm.cfg: Multiprocessing switched ON according to the parameter file. + 2020-03-25 12:31:52: oggm.cfg: Multiprocessing: using all available processors (N=8) + 2020-03-25 12:31:53: __main__: Starting OGGM inversion run + 2020-03-25 12:31:53: __main__: Number of glaciers: 35 + 2020-03-25 12:31:53: oggm.workflow: init_glacier_directories from prepro level 3 on 35 glaciers. + 2020-03-25 12:31:55: oggm.workflow: Execute entity task gdir_from_prepro on 35 glaciers + 2020-03-25 12:32:00: oggm.workflow: Execute entity task mass_conservation_inversion on 35 glaciers (...) - 2019-02-16 22:16:25: oggm.workflow: Execute entity task filter_inversion_output on 35 glaciers - 2019-02-16 22:16:25: oggm.workflow: Execute entity task glacier_statistics on 35 glaciers - 2019-02-16 22:16:25: __main__: OGGM is done! Time needed: 0:00:05 + 2020-03-25 12:32:33: oggm.workflow: Execute entity task filter_inversion_output on 35 glaciers + 2020-03-25 12:32:33: oggm.workflow: Execute entity task glacier_statistics on 35 glaciers + 2020-03-25 12:32:34: __main__: OGGM is done! Time needed: 0:00:41 Analysing the output diff --git a/docs/run_examples/run_rgi_region.rst b/docs/run_examples/run_rgi_region.rst index e9bb03c65..231f254c1 100644 --- a/docs/run_examples/run_rgi_region.rst +++ b/docs/run_examples/run_rgi_region.rst @@ -29,16 +29,19 @@ Script If everything went well, you should see an output similar to:: - 2019-02-16 17:50:51: oggm.cfg: Using configuration file: /home/mowglie/Documents/git/oggm-fork/oggm/params.cfg - 2019-02-16 17:50:52: __main__: Starting OGGM run - 2019-02-16 17:50:52: __main__: Number of glaciers: 54 - 2019-02-16 17:50:52: oggm.workflow: init_glacier_regions from prepro level 4 on 54 glaciers. - 2019-02-16 17:50:52: oggm.workflow: Execute entity task gdir_from_prepro on 54 glaciers - 2019-02-16 17:50:52: oggm.workflow: Multiprocessing: using all available processors (N=8) - 2019-02-16 17:50:54: oggm.workflow: Execute entity task run_random_climate on 54 glaciers - 2019-02-16 17:51:44: oggm.workflow: Execute entity task run_random_climate on 54 glaciers - 2019-02-16 17:52:36: oggm.workflow: Execute entity task run_random_climate on 54 glaciers - 2019-02-16 17:54:11: __main__: OGGM is done! Time needed: 0:03:20 + 2020-03-25 12:34:22: oggm.cfg: Using configuration file: /home/mowglie/disk/Dropbox/HomeDocs/git/oggm-fork/oggm/params.cfg + 2020-03-25 12:34:22: oggm.cfg: Multiprocessing switched ON according to the parameter file. + 2020-03-25 12:34:22: oggm.cfg: Multiprocessing: using all available processors (N=8) + 2020-03-25 12:34:23: __main__: Starting OGGM run + 2020-03-25 12:34:23: __main__: Number of glaciers: 54 + 2020-03-25 12:34:23: oggm.workflow: init_glacier_directories from prepro level 4 on 54 glaciers. + 2020-03-25 12:34:25: oggm.workflow: Execute entity task gdir_from_prepro on 54 glaciers + 2020-03-25 12:34:38: oggm.workflow: Execute entity task run_random_climate on 54 glaciers + 2020-03-25 12:35:22: oggm.workflow: Execute entity task run_random_climate on 54 glaciers + 2020-03-25 12:36:12: oggm.workflow: Execute entity task run_random_climate on 54 glaciers + 2020-03-25 12:38:03: oggm.workflow: Execute entity task glacier_statistics on 54 glaciers + 2020-03-25 12:38:04: __main__: OGGM is done! Time needed: 0:03:42 + Some analyses ------------- diff --git a/docs/whats-new.rst b/docs/whats-new.rst index 80e91c926..9a76f4a90 100644 --- a/docs/whats-new.rst +++ b/docs/whats-new.rst @@ -19,6 +19,15 @@ Breaking changes temperature standard deviation (:pull:`978`). The previous default was wrong and should not be the default. By `Fabien Maussion `_ +- Added a new "glacier directory initialization" global task: + `init_glacier_directories` (:pull:`983`, :issue:`965`). It replaces + `init_glacier_regions` and covers all its functionality, except that + it does *not* process the DEM data (this was a confusing "feature" of + `init_glacier_regions`). The old task `init_glacier_regions` is officially + deprecated but without warnings for now. Since it is a very widely used + task, we prefer to deprecate it in a slow cycle: first, change the + documentation, deprecate later. + By `Fabien Maussion `_ Enhancements ~~~~~~~~~~~~ diff --git a/oggm/cfg.py b/oggm/cfg.py index d3e698cc2..0d730d73f 100644 --- a/oggm/cfg.py +++ b/oggm/cfg.py @@ -633,8 +633,8 @@ def set_intersects_db(path_or_gdf=None): """Set the glacier intersection database for OGGM to use. It is now set automatically by the - :func:`oggm.workflow.init_glacier_regions` task, but setting it manually - can be useful for a slightly faster run initialization. + :func:`oggm.workflow.init_glacier_directories` task, but setting it + manually can be useful for a slightly faster run initialization. See :func:`oggm.utils.get_rgi_intersects_region_file` for how to obtain such data. diff --git a/oggm/cli/benchmark.py b/oggm/cli/benchmark.py index 37e72d486..e0e1f6273 100644 --- a/oggm/cli/benchmark.py +++ b/oggm/cli/benchmark.py @@ -127,8 +127,8 @@ def run_benchmark(rgi_version=None, rgi_reg=None, border=None, # Initialize working directories start = time.time() - gdirs = workflow.init_glacier_regions(rgidf, reset=True, force=True) - _add_time_to_df(odf, 'init_glacier_regions', time.time()-start) + gdirs = workflow.init_glacier_directories(rgidf, reset=True, force=True) + _add_time_to_df(odf, 'init_glacier_directories', time.time()-start) # Pre-download other files just in case if test_crudir is None: @@ -139,6 +139,7 @@ def run_benchmark(rgi_version=None, rgi_reg=None, border=None, # Tasks task_list = [ + tasks.define_glacier_region, tasks.process_cru_data, tasks.glacier_masks, tasks.compute_centerlines, diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index d59727e18..204525873 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -247,7 +247,8 @@ def _time_log(): rgidf['DEM_SOURCE'] = dem_source.upper() # L1 - go - gdirs = workflow.init_glacier_regions(rgidf, reset=True, force=True) + gdirs = workflow.init_glacier_directories(rgidf, reset=True, force=True) + workflow.execute_entity_task(tasks.define_glacier_region, gdirs) # Glacier stats sum_dir = os.path.join(base_dir, 'L1', 'summary') diff --git a/oggm/core/gis.py b/oggm/core/gis.py index c03609af3..b891a9063 100644 --- a/oggm/core/gis.py +++ b/oggm/core/gis.py @@ -118,56 +118,6 @@ def gaussian_blur(in_array, size): return scipy.signal.fftconvolve(padded_array, g, mode='valid') -def multi_to_poly(geometry, gdir=None): - """Sometimes an RGI geometry is a multipolygon: this should not happen. - - Parameters - ---------- - geometry : shpg.Polygon or shpg.MultiPolygon - the geometry to check - gdir : GlacierDirectory, optional - for logging - - Returns - ------- - the corrected geometry - """ - - # Log - rid = gdir.rgi_id + ': ' if gdir is not None else '' - - if 'Multi' in geometry.type: - parts = np.array(geometry) - for p in parts: - assert p.type == 'Polygon' - areas = np.array([p.area for p in parts]) - parts = parts[np.argsort(areas)][::-1] - areas = areas[np.argsort(areas)][::-1] - - # First case (was RGIV4): - # let's assume that one poly is exterior and that - # the other polygons are in fact interiors - exterior = parts[0].exterior - interiors = [] - was_interior = 0 - for p in parts[1:]: - if parts[0].contains(p): - interiors.append(p.exterior) - was_interior += 1 - if was_interior > 0: - # We are done here, good - geometry = shpg.Polygon(exterior, interiors) - else: - # This happens for bad geometries. We keep the largest - geometry = parts[0] - if np.any(areas[1:] > (areas[0] / 4)): - log.warning('Geometry {} lost quite a chunk.'.format(rid)) - - if geometry.type != 'Polygon': - raise InvalidGeometryError('Geometry {} is not a Polygon.'.format(rid)) - return geometry - - def _interp_polygon(polygon, dx): """Interpolates an irregular polygon to a regular step dx. @@ -280,70 +230,22 @@ def define_glacier_region(gdir, entity=None): gdir : :py:class:`oggm.GlacierDirectory` where to write the data entity : geopandas.GeoSeries - the glacier geometry to process + the glacier geometry to process - DEPRECATED. It is now ignored """ - # Make a local glacier map - proj_params = dict(name='tmerc', lat_0=0., lon_0=gdir.cenlon, - k=0.9996, x_0=0, y_0=0, datum='WGS84') - proj4_str = "+proj={name} +lat_0={lat_0} +lon_0={lon_0} +k={k} " \ - "+x_0={x_0} +y_0={y_0} +datum={datum}".format(**proj_params) - proj_in = pyproj.Proj("epsg:4326", preserve_units=True) - proj_out = pyproj.Proj(proj4_str, preserve_units=True) - project = partial(transform_proj, proj_in, proj_out) - # transform geometry to map - geometry = shapely.ops.transform(project, entity['geometry']) - geometry = multi_to_poly(geometry, gdir=gdir) - xx, yy = geometry.exterior.xy - - # Save transformed geometry to disk - entity = entity.copy() - entity['geometry'] = geometry - # Avoid fiona bug: https://github.com/Toblerity/Fiona/issues/365 - for k, s in entity.iteritems(): - if type(s) in [np.int32, np.int64]: - entity[k] = int(s) - towrite = gpd.GeoDataFrame(entity).T - towrite.crs = proj4_str - # Delete the source before writing - if 'DEM_SOURCE' in towrite: - del towrite['DEM_SOURCE'] + # Get the local map proj params and glacier extent + gdf = gdir.read_shapefile('outlines') + + # Get the map proj + utm_proj = salem.check_crs(gdf.crs) + + # Get glacier extent + xx, yy = gdf.iloc[0]['geometry'].exterior.xy # Define glacier area to use - area = entity['Area'] - - # Do we want to use the RGI area or ours? - if not cfg.PARAMS['use_rgi_area']: - # Update Area - area = geometry.area * 1e-6 - entity['Area'] = area - towrite['Area'] = area - - # Write shapefile - gdir.write_shapefile(towrite, 'outlines') - - # Also transform the intersects if necessary - gdf = cfg.PARAMS['intersects_gdf'] - if len(gdf) > 0: - gdf = gdf.loc[((gdf.RGIId_1 == gdir.rgi_id) | - (gdf.RGIId_2 == gdir.rgi_id))] - if len(gdf) > 0: - gdf = salem.transform_geopandas(gdf, to_crs=proj_out) - if hasattr(gdf.crs, 'srs'): - # salem uses pyproj - gdf.crs = gdf.crs.srs - gdir.write_shapefile(gdf, 'intersects') - else: - # Sanity check - if cfg.PARAMS['use_intersects']: - raise InvalidParamsError('You seem to have forgotten to set the ' - 'intersects file for this run. OGGM ' - 'works better with such a file. If you ' - 'know what your are doing, set ' - "cfg.PARAMS['use_intersects'] = False to " - "suppress this error.") - - # 6. choose a spatial resolution with respect to the glacier area + area = gdir.rgi_area_km2 + + # Choose a spatial resolution with respect to the glacier area dxmethod = cfg.PARAMS['grid_dx_method'] if dxmethod == 'linear': dx = np.rint(cfg.PARAMS['d1'] * area + cfg.PARAMS['d2']) @@ -382,12 +284,12 @@ def define_glacier_region(gdir, entity=None): ny = np.int((uly - lry) / dx) # Back to lon, lat for DEM download/preparation - tmp_grid = salem.Grid(proj=proj_out, nxny=(nx, ny), x0y0=(ulx, uly), + tmp_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), x0y0=(ulx, uly), dxdy=(dx, -dx), pixel_ref='corner') - minlon, maxlon, minlat, maxlat = tmp_grid.extent_in_crs(crs=salem.wgs84) # Open DEM + entity = gdf.iloc[0] source = entity.DEM_SOURCE if hasattr(entity, 'DEM_SOURCE') else None # We test DEM availability for glacier only (maps can grow big) if not is_dem_source_available(source, *gdir.extent_ll): @@ -431,7 +333,7 @@ def _get_nodata(rio_ds): # Set up profile for writing output profile = dem_dss[0].profile profile.update({ - 'crs': proj4_str, + 'crs': utm_proj.srs, 'transform': dst_transform, 'nodata': nodata, 'width': nx, @@ -463,7 +365,7 @@ def _get_nodata(rio_ds): # Destination parameters destination=dst_array, dst_transform=dst_transform, - dst_crs=proj4_str, + dst_crs=utm_proj.srs, dst_nodata=nodata, # Configuration resampling=resampling) @@ -474,7 +376,7 @@ def _get_nodata(rio_ds): # Glacier grid x0y0 = (ulx+dx/2, uly-dx/2) # To pixel center coordinates - glacier_grid = salem.Grid(proj=proj_out, nxny=(nx, ny), dxdy=(dx, -dx), + glacier_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), dxdy=(dx, -dx), x0y0=x0y0) glacier_grid.to_json(gdir.get_filepath('glacier_grid')) diff --git a/oggm/tests/funcs.py b/oggm/tests/funcs.py index d1479d1c8..b2495cba2 100644 --- a/oggm/tests/funcs.py +++ b/oggm/tests/funcs.py @@ -310,7 +310,7 @@ def init_hef(reset=False, border=40, logging_level='INFO'): if not reset: return gdir - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) execute_entity_task(gis.glacier_masks, [gdir]) execute_entity_task(centerlines.compute_centerlines, [gdir]) centerlines.initialize_flowlines(gdir) @@ -390,7 +390,7 @@ def init_columbia(reset=False): if gdir.has_file('climate_monthly'): return gdir - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) diff --git a/oggm/tests/test_benchmarks.py b/oggm/tests/test_benchmarks.py index e2b7a3a62..098a5f2e8 100644 --- a/oggm/tests/test_benchmarks.py +++ b/oggm/tests/test_benchmarks.py @@ -97,10 +97,11 @@ def test_mb(self): rgidf = gpd.read_file(get_demo_file('SouthGlacier.shp')) # Go - initialize working directories - gdirs = workflow.init_glacier_regions(rgidf) + gdirs = workflow.init_glacier_directories(rgidf) # Preprocessing tasks task_list = [ + tasks.define_glacier_region, tasks.glacier_masks, tasks.compute_centerlines, tasks.initialize_flowlines, @@ -159,10 +160,11 @@ def test_inversion_attributes(self): rgidf = gpd.read_file(get_demo_file('SouthGlacier.shp')) # Go - initialize working directories - gdirs = workflow.init_glacier_regions(rgidf) + gdirs = workflow.init_glacier_directories(rgidf) # Preprocessing tasks task_list = [ + tasks.define_glacier_region, tasks.glacier_masks, tasks.compute_centerlines, tasks.initialize_flowlines, @@ -237,10 +239,11 @@ def test_inversion(self): rgidf = gpd.read_file(get_demo_file('SouthGlacier.shp')) # Go - initialize working directories - gdirs = workflow.init_glacier_regions(rgidf) + gdirs = workflow.init_glacier_directories(rgidf) # Preprocessing tasks task_list = [ + tasks.define_glacier_region, tasks.glacier_masks, tasks.compute_centerlines, tasks.initialize_flowlines, @@ -307,10 +310,11 @@ def test_optimize_inversion(self): rgidf = gpd.read_file(get_demo_file('SouthGlacier.shp')) # Go - initialize working directories - gdirs = workflow.init_glacier_regions(rgidf) + gdirs = workflow.init_glacier_directories(rgidf) # Preprocessing tasks task_list = [ + tasks.define_glacier_region, tasks.glacier_masks, tasks.compute_centerlines, tasks.initialize_flowlines, @@ -385,10 +389,11 @@ def test_workflow(self): rgidf = gpd.read_file(get_demo_file('SouthGlacier.shp')) # Go - initialize working directories - gdirs = workflow.init_glacier_regions(rgidf) + gdirs = workflow.init_glacier_directories(rgidf) # Preprocessing tasks task_list = [ + tasks.define_glacier_region, tasks.glacier_masks, tasks.compute_centerlines, tasks.initialize_flowlines, @@ -453,7 +458,7 @@ def test_set_width(self): entity = gpd.read_file(self.rgi_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -510,7 +515,7 @@ def test_run(self): entity = gpd.read_file(self.rgi_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) diff --git a/oggm/tests/test_graphics.py b/oggm/tests/test_graphics.py index 054eba3ec..0bfd4f833 100644 --- a/oggm/tests/test_graphics.py +++ b/oggm/tests/test_graphics.py @@ -163,7 +163,7 @@ def test_multiple_inversion(): hef_rgi = gpd.read_file(get_demo_file('divides_hef.shp')) hef_rgi.loc[0, 'RGIId'] = 'RGI50-11.00897' - gdirs = workflow.init_glacier_regions(hef_rgi) + gdirs = workflow.init_glacier_directories(hef_rgi) workflow.gis_prepro_tasks(gdirs) workflow.climate_tasks(gdirs) workflow.inversion_tasks(gdirs) @@ -241,7 +241,7 @@ def test_multiple_models(): hef_rgi = gpd.read_file(get_demo_file('divides_hef.shp')) hef_rgi.loc[0, 'RGIId'] = 'RGI50-11.00897' - gdirs = workflow.init_glacier_regions(hef_rgi) + gdirs = workflow.init_glacier_directories(hef_rgi) workflow.gis_prepro_tasks(gdirs) workflow.climate_tasks(gdirs) workflow.inversion_tasks(gdirs) @@ -311,7 +311,7 @@ def test_chhota_shigri(): df['Area'] = df.Area * 1e-6 # cause it was in m2 df['RGIId'] = ['RGI50-14.15990' + d for d in ['_d01', '_d02']] - gdirs = workflow.init_glacier_regions(df) + gdirs = workflow.init_glacier_directories(df) workflow.gis_prepro_tasks(gdirs) for gdir in gdirs: climate.apparent_mb_from_linear_mb(gdir) @@ -350,7 +350,7 @@ def test_ice_cap(): df['Area'] = df.Area * 1e-6 # cause it was in m2 df['RGIId'] = ['RGI50-05.08389_d{:02d}'.format(d+1) for d in df.index] - gdirs = workflow.init_glacier_regions(df) + gdirs = workflow.init_glacier_directories(df) workflow.gis_prepro_tasks(gdirs) from salem import mercator_grid, Map @@ -386,7 +386,7 @@ def test_coxe(): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=testdir, reset=True) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) diff --git a/oggm/tests/test_models.py b/oggm/tests/test_models.py index 11d5febbb..4af1b0d90 100644 --- a/oggm/tests/test_models.py +++ b/oggm/tests/test_models.py @@ -182,7 +182,7 @@ def test_define_divides(self, class_case_dir): # This is another glacier with divides entity = rgidf.loc[rgidf.RGIId == 'RGI50-11.00719_d01'].iloc[0] gdir = GlacierDirectory(entity, base_dir=class_case_dir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -1597,7 +1597,7 @@ def inversion_gdir(class_case_dir): entity = gpd.read_file(hef_file).iloc[0] gdir = GlacierDirectory(entity, base_dir=class_case_dir, reset=True) - define_glacier_region(gdir, entity=entity) + define_glacier_region(gdir) return gdir @@ -2844,7 +2844,7 @@ def test_merged_simulation(self): (rgidf.RGIId == 'RGI50-11.00719_d01') | (rgidf.RGIId == 'RGI50-11.00779') | (rgidf.RGIId == 'RGI50-11.00746')].copy() - gdirs = workflow.init_glacier_regions(glcdf) + gdirs = workflow.init_glacier_directories(glcdf) workflow.gis_prepro_tasks(gdirs) workflow.climate_tasks(gdirs) workflow.inversion_tasks(gdirs) diff --git a/oggm/tests/test_prepro.py b/oggm/tests/test_prepro.py index 3552c233d..59bb78572 100644 --- a/oggm/tests/test_prepro.py +++ b/oggm/tests/test_prepro.py @@ -94,13 +94,15 @@ def test_define_region(self): hef_file = get_demo_file('Hintereisferner_RGI5.shp') entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) extent = gdir.extent_ll tdf = gdir.read_shapefile('outlines') myarea = tdf.geometry.area * 10**-6 np.testing.assert_allclose(myarea, np.float(tdf['Area']), rtol=1e-2) self.assertTrue(gdir.has_file('intersects')) + np.testing.assert_array_equal(gdir.intersects_ids, + ['RGI50-11.00846', 'RGI50-11.00950']) # From string gdir = oggm.GlacierDirectory(gdir.rgi_id, base_dir=self.testdir) @@ -115,7 +117,7 @@ def test_define_region(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir, reset=True) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) # Close but not same assert gdir.rgi_area_km2 != prev_area assert gdir.cenlon != prev_lon @@ -132,8 +134,7 @@ def test_init_glacier_regions(self): gdir = workflow.init_glacier_regions(hef_rgi)[0] nx, ny = gdir.grid.nx, gdir.grid.ny - # Change something and note that no error happens because the - # directory is not overwritten + # Change something and note that no change occurs because dem is there cfg.PARAMS['border'] = 12 gdir = workflow.init_glacier_regions(hef_rgi)[0] assert nx == gdir.grid.nx @@ -149,7 +150,7 @@ def test_divides_as_glaciers(self): ['_d01', '_d02', '_d03']] # Just check that things are working - gdirs = workflow.init_glacier_regions(hef_rgi) + gdirs = workflow.init_glacier_directories(hef_rgi) workflow.gis_prepro_tasks(gdirs) assert gdirs[0].rgi_id == 'RGI50-11.00897_d01' @@ -163,7 +164,7 @@ def test_dx_methods(self): # Test fixed method cfg.PARAMS['grid_dx_method'] = 'fixed' cfg.PARAMS['fixed_dx'] = 50 - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) mygrid = salem.Grid.from_json(gdir.get_filepath('glacier_grid')) np.testing.assert_allclose(np.abs(mygrid.dx), 50.) @@ -172,7 +173,7 @@ def test_dx_methods(self): cfg.PARAMS['d1'] = 5. cfg.PARAMS['d2'] = 10. cfg.PARAMS['dmax'] = 100. - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) targetdx = np.rint(5. * gdir.rgi_area_km2 + 10.) targetdx = np.clip(targetdx, 10., 100.) mygrid = salem.Grid.from_json(gdir.get_filepath('glacier_grid')) @@ -183,7 +184,7 @@ def test_dx_methods(self): cfg.PARAMS['d1'] = 5. cfg.PARAMS['d2'] = 10. cfg.PARAMS['dmax'] = 100. - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) targetdx = np.rint(5. * np.sqrt(gdir.rgi_area_km2) + 10.) targetdx = np.clip(targetdx, 10., 100.) mygrid = salem.Grid.from_json(gdir.get_filepath('glacier_grid')) @@ -209,7 +210,7 @@ def test_repr(self): hef_file = get_demo_file('Hintereisferner_RGI5.shp') entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) self.assertEqual(gdir.__repr__(), expected) def test_glacierdir(self): @@ -217,7 +218,7 @@ def test_glacierdir(self): hef_file = get_demo_file('Hintereisferner_RGI5.shp') entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) # this should simply run oggm.GlacierDirectory(entity.RGIId, base_dir=self.testdir) @@ -229,7 +230,7 @@ def test_glacier_masks(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.process_dem(gdir) gis.glacier_masks(gdir) gis.gridded_attributes(gdir) @@ -265,7 +266,7 @@ def test_simple_glacier_masks(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.simple_glacier_masks(gdir) with utils.ncDataset(gdir.get_filepath('gridded_data')) as nc: @@ -324,7 +325,7 @@ def test_glacier_masks_other_glacier(self): cfg.PARAMS['border'] = 1 gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) # The test below does NOT pass on OGGM shutil.copyfile(gdir.get_filepath('gridded_data'), @@ -354,7 +355,7 @@ def test_rasterio_glacier_masks(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) # specifying a source will look for a DEN in a respective folder self.assertRaises(ValueError, gis.rasterio_glacier_mask, @@ -394,7 +395,7 @@ def test_intersects(self): hef_file = get_demo_file('Hintereisferner_RGI5.shp') entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) self.assertTrue(gdir.has_file('intersects')) def test_dem_source_text(self): @@ -407,7 +408,7 @@ def test_dem_daterange_dateinfo(self): hef_file = get_demo_file('Hintereisferner_RGI5.shp') entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) # dem_info should return a string self.assertIsInstance(gdir.dem_info, str) @@ -431,7 +432,7 @@ def test_custom_basename(self): hef_file = get_demo_file('Hintereisferner_RGI5.shp') entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) cfg.add_to_basenames('mybn', 'testfb.pkl', docstr='Some docs') @@ -528,7 +529,7 @@ def test_centerlines(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) @@ -552,7 +553,7 @@ def test_downstream(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -573,7 +574,7 @@ def test_downstream_bedshape(self): cfg.PARAMS['border'] = 80 gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -628,26 +629,29 @@ def test_baltoro_centerlines(self): cfg.PATHS['dem_file'] = get_demo_file('baltoro_srtm_clip.tif') b_file = get_demo_file('baltoro_wgs84.shp') - entity = gpd.read_file(b_file).iloc[0] + gdf = gpd.read_file(b_file) kienholz_file = get_demo_file('centerlines_baltoro_wgs84.shp') kdf = gpd.read_file(kienholz_file) # add fake attribs - del entity['RGIID'] - entity.RGIId = 'RGI50-00.00000' - entity['GLIMSId'] = entity['GLIMSID'] - entity['Area'] = entity['AREA'] - entity['CenLat'] = entity['CENLAT'] - entity['CenLon'] = entity['CENLON'] - entity.BgnDate = '-999' - entity.Name = 'Baltoro' - entity.GlacType = '0000' - entity.Status = '0' - entity.O1Region = '01' - entity.O2Region = '01' + area = gdf['AREA'] + del gdf['RGIID'] + del gdf['AREA'] + gdf['RGIId'] = 'RGI50-00.00000' + gdf['GLIMSId'] = gdf['GLIMSID'] + gdf['Area'] = area + gdf['CenLat'] = gdf['CENLAT'] + gdf['CenLon'] = gdf['CENLON'] + gdf['BgnDate'] = '-999' + gdf['Name'] = 'Baltoro' + gdf['GlacType'] = '0000' + gdf['Status'] = '0' + gdf['O1Region'] = '01' + gdf['O2Region'] = '01' + entity = gdf.iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) @@ -747,7 +751,7 @@ def test_catchment_area(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.catchment_area(gdir) @@ -772,7 +776,7 @@ def test_flowlines(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -807,7 +811,7 @@ def test_geom_width(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -821,7 +825,7 @@ def test_width(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -883,7 +887,7 @@ def test_nodivides_correct_slope(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -938,7 +942,7 @@ def test_distribute_climate(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) climate.process_custom_climate_data(gdir) ci = gdir.read_json('climate_info') @@ -963,7 +967,7 @@ def test_distribute_climate_grad(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) climate.process_custom_climate_data(gdir) ci = gdir.read_json('climate_info') @@ -984,7 +988,7 @@ def test_distribute_climate_parallel(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) climate.process_custom_climate_data(gdir) ci = gdir.read_json('climate_info') @@ -1010,10 +1014,10 @@ def test_distribute_climate_cru(self): gdirs = [] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gdirs.append(gdir) gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir_cru) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gdirs.append(gdir) climate.process_custom_climate_data(gdirs[0]) @@ -1053,10 +1057,10 @@ def test_distribute_climate_dummy(self): gdirs = [] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gdirs.append(gdir) gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir_cru) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gdirs.append(gdir) climate.process_dummy_cru_file(gdirs[0], seed=0) @@ -1102,10 +1106,10 @@ def test_distribute_climate_histalp_new(self): gdirs = [] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gdirs.append(gdir) gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir_cru) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gdirs.append(gdir) climate.process_custom_climate_data(gdirs[0]) @@ -1159,12 +1163,12 @@ def test_sh(self): assert gdir.hemisphere == 'nh' gdir.hemisphere = 'sh' - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gdirs.append(gdir) gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir_cru) assert gdir.hemisphere == 'nh' gdir.hemisphere = 'sh' - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gdirs.append(gdir) climate.process_custom_climate_data(gdirs[0]) ci = gdirs[0].read_json('climate_info') @@ -1217,7 +1221,7 @@ def test_mb_climate(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) climate.process_custom_climate_data(gdir) with utils.ncDataset(get_demo_file('histalp_merged_hef.nc')) as nc_r: @@ -1280,7 +1284,7 @@ def test_yearly_mb_climate(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) climate.process_custom_climate_data(gdir) with utils.ncDataset(get_demo_file('histalp_merged_hef.nc')) as nc_r: @@ -1360,7 +1364,7 @@ def test_mu_candidates(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -1393,7 +1397,7 @@ def test_find_tstars(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -1453,7 +1457,7 @@ def test_find_tstars_multiple_mus(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -1499,7 +1503,7 @@ def test_local_t_star(self): cfg.PARAMS['prcp_scaling_factor'] = 2.9 gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -1572,7 +1576,7 @@ def test_local_t_star_fallback(self): cfg.PARAMS['continue_on_error'] = True gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -1662,14 +1666,14 @@ def test_automated_workflow(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) assert gdir.rgi_version == '50' - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) workflow.gis_prepro_tasks([gdir]) workflow.climate_tasks([gdir]) hef_file = get_demo_file('Hintereisferner_RGI6.shp') entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) assert gdir.rgi_version == '60' workflow.gis_prepro_tasks([gdir]) workflow.climate_tasks([gdir]) @@ -1713,7 +1717,7 @@ def test_filter(self): cfg.PARAMS['correct_for_neg_flux'] = False gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -1742,7 +1746,7 @@ def test_correct(self): cfg.PARAMS['correct_for_neg_flux'] = False gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -1803,7 +1807,7 @@ def test_and_compare_two_methods(self): entity = entity.loc[entity.RGIId == 'RGI50-11.00666'].iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -1889,7 +1893,7 @@ def test_invert_hef(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -2015,7 +2019,7 @@ def test_invert_hef_from_linear_mb(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -2110,7 +2114,7 @@ def test_invert_hef_from_any_mb(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -2142,7 +2146,7 @@ def test_distribute(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -2203,7 +2207,7 @@ def test_invert_hef_nofs(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -2310,7 +2314,7 @@ def test_continue_on_error(self): entity.RGIId = 'RGI50-11.fake' gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -2364,7 +2368,7 @@ def test_inversion_and_run_with_calving(self): entity = gpd.read_file(coxe_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -2769,7 +2773,7 @@ def test_process_cesm(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) climate.process_cru_data(gdir) ci = gdir.read_json('climate_info') @@ -2824,7 +2828,7 @@ def test_process_cmip5(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) climate.process_cru_data(gdir) ci = gdir.read_json('climate_info') @@ -2882,7 +2886,7 @@ def test_process_cmip5_scale(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) climate.process_cru_data(gdir) ci = gdir.read_json('climate_info') @@ -2954,7 +2958,7 @@ def test_compile_climate_input(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) climate.process_cru_data(gdir) utils.compile_climate_input([gdir]) @@ -3044,7 +3048,7 @@ def test_invert(self): entity = gpd.read_file(hef_file).iloc[0] gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) centerlines.compute_centerlines(gdir) centerlines.initialize_flowlines(gdir) @@ -3113,7 +3117,7 @@ def test_pipe_log(self): cfg.PARAMS['continue_on_error'] = True gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) # This will "run" but log an error @@ -3137,7 +3141,7 @@ def test_task_status(self): cfg.PARAMS['continue_on_error'] = True gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) gis.glacier_masks(gdir) self.assertEqual(gdir.get_task_status(gis.glacier_masks.__name__), @@ -3239,7 +3243,7 @@ def test_read_gmip_data(self): entity = gpd.read_file(hef_file).iloc[0] entity['RGIId'] = 'RGI60-11.00897' gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) from oggm.sandbox import pygem_compat area_path = get_demo_file('gmip_area_centraleurope_10_sel.dat') @@ -3258,7 +3262,7 @@ def test_flowlines_from_gmip_data(self): entity = gpd.read_file(hef_file).iloc[0] entity['RGIId'] = 'RGI60-11.00897' gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) - gis.define_glacier_region(gdir, entity=entity) + gis.define_glacier_region(gdir) from oggm.sandbox import pygem_compat area_path = get_demo_file('gmip_area_centraleurope_10_sel.dat') diff --git a/oggm/tests/test_utils.py b/oggm/tests/test_utils.py index 429188880..195adc033 100644 --- a/oggm/tests/test_utils.py +++ b/oggm/tests/test_utils.py @@ -355,6 +355,7 @@ def test_leclercq_data(self): hef_file = utils.get_demo_file('Hintereisferner_RGI5.shp') entity = gpd.read_file(hef_file).iloc[0] + cfg.PARAMS['use_intersects'] = False gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir) df = gdir.get_ref_length_data() @@ -421,14 +422,15 @@ def clean_dir(self): def test_to_and_from_tar(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(self.rgidf) + gdirs = workflow.init_glacier_directories(self.rgidf) + workflow.execute_entity_task(tasks.define_glacier_region, gdirs) # End - compress all workflow.execute_entity_task(utils.gdir_to_tar, gdirs) # Test - reopen form tar - gdirs = workflow.init_glacier_regions(self.rgidf, from_tar=True, - delete_tar=True) + gdirs = workflow.init_glacier_directories(self.rgidf, from_tar=True, + delete_tar=True) for gdir in gdirs: assert gdir.has_file('dem') assert not os.path.exists(gdir.dir + '.tar.gz') @@ -436,7 +438,7 @@ def test_to_and_from_tar(self): workflow.execute_entity_task(utils.gdir_to_tar, gdirs) - gdirs = workflow.init_glacier_regions(self.rgidf, from_tar=True) + gdirs = workflow.init_glacier_directories(self.rgidf, from_tar=True) for gdir in gdirs: assert gdir.has_file('gridded_data') assert os.path.exists(gdir.dir + '.tar.gz') @@ -444,14 +446,15 @@ def test_to_and_from_tar(self): def test_to_and_from_basedir_tar(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(self.rgidf) + gdirs = workflow.init_glacier_directories(self.rgidf) + workflow.execute_entity_task(tasks.define_glacier_region, gdirs) # End - compress all workflow.execute_entity_task(utils.gdir_to_tar, gdirs) utils.base_dir_to_tar() # Test - reopen form tar - gdirs = workflow.init_glacier_regions(self.rgidf, from_tar=True) + gdirs = workflow.init_glacier_directories(self.rgidf, from_tar=True) for gdir in gdirs: assert gdir.has_file('dem') @@ -460,7 +463,7 @@ def test_to_and_from_basedir_tar(self): workflow.execute_entity_task(utils.gdir_to_tar, gdirs) - gdirs = workflow.init_glacier_regions(self.rgidf, from_tar=True) + gdirs = workflow.init_glacier_directories(self.rgidf, from_tar=True) for gdir in gdirs: assert gdir.has_file('gridded_data') assert os.path.exists(gdir.dir + '.tar.gz') @@ -468,7 +471,8 @@ def test_to_and_from_basedir_tar(self): def test_to_and_from_tar_new_dir(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(self.rgidf) + gdirs = workflow.init_glacier_directories(self.rgidf) + workflow.execute_entity_task(tasks.define_glacier_region, gdirs) # End - compress all base_dir = os.path.join(self.testdir, 'new_base_dir') @@ -480,8 +484,9 @@ def test_to_and_from_tar_new_dir(self): assert base_dir in p shutil.copyfile(p, os.path.normpath(gdir.dir) + '.tar.gz') - gdirs = workflow.init_glacier_regions(self.rgidf, from_tar=True, - delete_tar=True) + gdirs = workflow.init_glacier_directories(self.rgidf, + from_tar=True, + delete_tar=True) for gdir in gdirs: assert gdir.has_file('dem') assert not os.path.exists(gdir.dir + '.tar.gz') @@ -491,7 +496,8 @@ def test_to_and_from_tar_new_dir(self): def test_to_and_from_tar_string(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(self.rgidf) + gdirs = workflow.init_glacier_directories(self.rgidf) + workflow.execute_entity_task(tasks.define_glacier_region, gdirs) # End - compress all base_dir = os.path.join(self.testdir, 'new_base_dir') @@ -546,11 +552,11 @@ def clean_dir(self): def test_start_from_level_1(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(self.rgidf.iloc[:2], - from_prepro_level=1, - prepro_rgi_version='61', - prepro_border=20, - use_demo_glaciers=False) + gdirs = workflow.init_glacier_directories(self.rgidf.iloc[:2], + from_prepro_level=1, + prepro_rgi_version='61', + prepro_border=20, + use_demo_glaciers=False) n_intersects = 0 for gdir in gdirs: assert gdir.has_file('dem') @@ -564,9 +570,9 @@ def test_start_from_level_1_str(self): # Go - initialize working directories entitites = self.rgidf.iloc[:2].RGIId cfg.PARAMS['border'] = 20 - gdirs = workflow.init_glacier_regions(entitites, - from_prepro_level=1, - use_demo_glaciers=False) + gdirs = workflow.init_glacier_directories(entitites, + from_prepro_level=1, + use_demo_glaciers=False) n_intersects = 0 for gdir in gdirs: assert gdir.has_file('dem') @@ -576,9 +582,9 @@ def test_start_from_level_1_str(self): # One string cfg.PARAMS['border'] = 20 - gdirs = workflow.init_glacier_regions('RGI60-11.00897', - from_prepro_level=1, - use_demo_glaciers=False) + gdirs = workflow.init_glacier_directories('RGI60-11.00897', + from_prepro_level=1, + use_demo_glaciers=False) n_intersects = 0 for gdir in gdirs: assert gdir.has_file('dem') @@ -590,11 +596,11 @@ def test_start_from_level_1_str(self): def test_start_from_level_2(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(self.rgidf.iloc[:2], - from_prepro_level=2, - prepro_rgi_version='61', - prepro_border=20, - use_demo_glaciers=False) + gdirs = workflow.init_glacier_directories(self.rgidf.iloc[:2], + from_prepro_level=2, + prepro_rgi_version='61', + prepro_border=20, + use_demo_glaciers=False) n_intersects = 0 for gdir in gdirs: assert gdir.has_file('dem') @@ -607,11 +613,11 @@ def test_start_from_level_2(self): def test_start_from_level_3(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(self.rgidf.iloc[:2], - from_prepro_level=3, - prepro_rgi_version='61', - prepro_border=20, - use_demo_glaciers=False) + gdirs = workflow.init_glacier_directories(self.rgidf.iloc[:2], + from_prepro_level=3, + prepro_rgi_version='61', + prepro_border=20, + use_demo_glaciers=False) n_intersects = 0 for gdir in gdirs: assert gdir.has_file('dem') @@ -645,21 +651,21 @@ def test_start_from_level_3(self): def test_start_from_level_4(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(self.rgidf.iloc[:2], - from_prepro_level=4, - prepro_rgi_version='61', - prepro_border=20, - use_demo_glaciers=False) + gdirs = workflow.init_glacier_directories(self.rgidf.iloc[:2], + from_prepro_level=4, + prepro_rgi_version='61', + prepro_border=20, + use_demo_glaciers=False) workflow.execute_entity_task(tasks.run_random_climate, gdirs, nyears=10) def test_start_from_demo(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(['kwf', 'hef'], - from_prepro_level=4, - prepro_rgi_version='61', - prepro_border=10) + gdirs = workflow.init_glacier_directories(['kwf', 'hef'], + from_prepro_level=4, + prepro_rgi_version='61', + prepro_border=10) workflow.execute_entity_task(tasks.run_random_climate, gdirs, nyears=10) @@ -667,10 +673,10 @@ def test_start_from_demo(self): def test_corrupted_file(self): # Go - initialize working directories - gdirs = workflow.init_glacier_regions(['hef'], - from_prepro_level=4, - prepro_rgi_version='61', - prepro_border=10) + gdirs = workflow.init_glacier_directories(['hef'], + from_prepro_level=4, + prepro_rgi_version='61', + prepro_border=10) cfile = utils.get_prepro_gdir('61', 'RGI60-11.00787', 10, 4, base_url=utils.DEMO_GDIR_URL) @@ -682,10 +688,10 @@ def test_corrupted_file(self): f.write('ups') # This should retrigger a download and just work - gdirs = workflow.init_glacier_regions(['hef'], - from_prepro_level=4, - prepro_rgi_version='61', - prepro_border=10) + gdirs = workflow.init_glacier_directories(['hef'], + from_prepro_level=4, + prepro_rgi_version='61', + prepro_border=10) workflow.execute_entity_task(tasks.run_random_climate, gdirs, nyears=10) @@ -2358,11 +2364,11 @@ def test_from_prepro(self): for rid in rgidf.RGIId] # Go - initialize working directories - gdirs = workflow.init_glacier_regions(rgidf.iloc[:2], - from_prepro_level=1, - prepro_rgi_version='61', - prepro_border=10, - use_demo_glaciers=False) + gdirs = workflow.init_glacier_directories(rgidf.iloc[:2], + from_prepro_level=1, + prepro_rgi_version='61', + prepro_border=10, + use_demo_glaciers=False) n_intersects = 0 for gdir in gdirs: assert gdir.has_file('dem') diff --git a/oggm/tests/test_workflow.py b/oggm/tests/test_workflow.py index f5f9c0b8d..a66585961 100644 --- a/oggm/tests/test_workflow.py +++ b/oggm/tests/test_workflow.py @@ -88,7 +88,7 @@ def up_to_climate(reset=False): cfg.PARAMS['run_mb_calibration'] = True # Go - gdirs = workflow.init_glacier_regions(rgidf) + gdirs = workflow.init_glacier_directories(rgidf) try: tasks.catchment_width_correction(gdirs[0]) diff --git a/oggm/utils/_funcs.py b/oggm/utils/_funcs.py index dff196eee..65c1d0b2f 100644 --- a/oggm/utils/_funcs.py +++ b/oggm/utils/_funcs.py @@ -28,7 +28,7 @@ import oggm.cfg as cfg from oggm.cfg import SEC_IN_YEAR, SEC_IN_MONTH from oggm.utils._downloads import get_demo_file -from oggm.exceptions import InvalidParamsError +from oggm.exceptions import InvalidParamsError, InvalidGeometryError # Module logger log = logging.getLogger('.'.join(__name__.split('.')[:-1])) @@ -467,6 +467,56 @@ def polygon_intersections(gdf): return out +def multipolygon_to_polygon(geometry, gdir=None): + """Sometimes an RGI geometry is a multipolygon: this should not happen. + + Parameters + ---------- + geometry : shpg.Polygon or shpg.MultiPolygon + the geometry to check + gdir : GlacierDirectory, optional + for logging + + Returns + ------- + the corrected geometry + """ + + # Log + rid = gdir.rgi_id + ': ' if gdir is not None else '' + + if 'Multi' in geometry.type: + parts = np.array(geometry) + for p in parts: + assert p.type == 'Polygon' + areas = np.array([p.area for p in parts]) + parts = parts[np.argsort(areas)][::-1] + areas = areas[np.argsort(areas)][::-1] + + # First case (was RGIV4): + # let's assume that one poly is exterior and that + # the other polygons are in fact interiors + exterior = parts[0].exterior + interiors = [] + was_interior = 0 + for p in parts[1:]: + if parts[0].contains(p): + interiors.append(p.exterior) + was_interior += 1 + if was_interior > 0: + # We are done here, good + geometry = shpg.Polygon(exterior, interiors) + else: + # This happens for bad geometries. We keep the largest + geometry = parts[0] + if np.any(areas[1:] > (areas[0] / 4)): + log.warning('Geometry {} lost quite a chunk.'.format(rid)) + + if geometry.type != 'Polygon': + raise InvalidGeometryError('Geometry {} is not a Polygon.'.format(rid)) + return geometry + + def floatyear_to_date(yr): """Converts a float year to an actual (year, month) pair. diff --git a/oggm/utils/_workflow.py b/oggm/utils/_workflow.py index a4f86a488..0192c23c1 100644 --- a/oggm/utils/_workflow.py +++ b/oggm/utils/_workflow.py @@ -41,14 +41,20 @@ try: import salem from salem import wgs84 + from salem.gis import transform_proj +except ImportError: + pass +try: + import pyproj except ImportError: pass + # Locals from oggm import __version__ from oggm.utils._funcs import (calendardate_to_hydrodate, date_to_floatyear, tolist, filter_rgi_name, parse_rgi_meta, - haversine) + haversine, multipolygon_to_polygon) from oggm.utils._downloads import (get_demo_file, get_wgms_files, get_rgi_glacier_entities) from oggm import cfg @@ -1564,10 +1570,12 @@ def __init__(self, rgi_entity, base_dir=None, reset=False, xx, yy = salem.transform_proj(crs, salem.wgs84, [g.bounds[0], g.bounds[2]], [g.bounds[1], g.bounds[3]]) + write_shp = False else: g = rgi_entity['geometry'] xx, yy = ([g.bounds[0], g.bounds[2]], [g.bounds[1], g.bounds[3]]) + write_shp = True # Extent of the glacier in lon/lat self.extent_ll = [xx, yy] @@ -1680,6 +1688,7 @@ def __init__(self, rgi_entity, base_dir=None, reset=False, if from_tar is True: from_tar = self.dir + '.tar.gz' robust_tar_extract(from_tar, self.dir, delete_tar=delete_tar) + write_shp = False else: mkdir(self.dir) @@ -1689,6 +1698,10 @@ def __init__(self, rgi_entity, base_dir=None, reset=False, # logging file self.logfile = os.path.join(self.dir, 'log.txt') + if write_shp: + # Write shapefile + self._reproject_and_write_shapefile(rgi_entity) + # Optimization self._mbdf = None self._mbprofdf = None @@ -1713,6 +1726,68 @@ def __repr__(self): str(self.grid.dy) + ')'] return '\n'.join(summary) + '\n' + def _reproject_and_write_shapefile(self, entity): + + # Make a local glacier map + params = dict(name='tmerc', lat_0=0., lon_0=self.cenlon, + k=0.9996, x_0=0, y_0=0, datum='WGS84') + proj4_str = "+proj={name} +lat_0={lat_0} +lon_0={lon_0} +k={k} " \ + "+x_0={x_0} +y_0={y_0} +datum={datum}".format(**params) + + # Reproject + proj_in = pyproj.Proj("epsg:4326", preserve_units=True) + proj_out = pyproj.Proj(proj4_str, preserve_units=True) + project = partial(transform_proj, proj_in, proj_out) + + # transform geometry to map + geometry = shp_trafo(project, entity['geometry']) + geometry = multipolygon_to_polygon(geometry, gdir=self) + + # Save transformed geometry to disk + entity = entity.copy() + entity['geometry'] = geometry + + # Do we want to use the RGI area or ours? + if not cfg.PARAMS['use_rgi_area']: + # Update Area + area = geometry.area * 1e-6 + entity['Area'] = area + + # Avoid fiona bug: https://github.com/Toblerity/Fiona/issues/365 + for k, s in entity.iteritems(): + if type(s) in [np.int32, np.int64]: + entity[k] = int(s) + towrite = gpd.GeoDataFrame(entity).T + towrite.crs = proj4_str + # Delete the source before writing + if 'DEM_SOURCE' in towrite: + del towrite['DEM_SOURCE'] + + # Write shapefile + self.write_shapefile(towrite, 'outlines') + + # Also transform the intersects if necessary + gdf = cfg.PARAMS['intersects_gdf'] + if len(gdf) > 0: + gdf = gdf.loc[((gdf.RGIId_1 == self.rgi_id) | + (gdf.RGIId_2 == self.rgi_id))] + if len(gdf) > 0: + gdf = salem.transform_geopandas(gdf, to_crs=proj_out) + if hasattr(gdf.crs, 'srs'): + # salem uses pyproj + gdf.crs = gdf.crs.srs + self.write_shapefile(gdf, 'intersects') + else: + # Sanity check + if cfg.PARAMS['use_intersects']: + raise InvalidParamsError( + 'You seem to have forgotten to set the ' + 'intersects file for this run. OGGM ' + 'works better with such a file. If you ' + 'know what your are doing, set ' + "cfg.PARAMS['use_intersects'] = False to " + "suppress this error.") + @lazy_property def grid(self): """A ``salem.Grid`` handling the georeferencing of the local grid""" @@ -1725,8 +1800,19 @@ def rgi_area_km2(self): _area = self.read_shapefile('outlines')['Area'] return np.round(float(_area), decimals=3) except OSError: - raise RuntimeError('Please run `define_glacier_region` before ' - 'using this property.') + raise RuntimeError('No outlines available') + + @lazy_property + def intersects_ids(self): + """The glacier's intersects RGI ids.""" + try: + gdf = self.read_shapefile('intersects') + ids = np.append(gdf['RGIId_1'], gdf['RGIId_2']) + ids = list(np.unique(np.sort(ids))) + ids.remove(self.rgi_id) + return ids + except OSError: + return [] @lazy_property def dem_daterange(self): diff --git a/oggm/workflow.py b/oggm/workflow.py index 4bd8f434a..a40c8bfaa 100644 --- a/oggm/workflow.py +++ b/oggm/workflow.py @@ -53,8 +53,6 @@ def init_mp_pool(reset=False): global_lock = mp.Manager().Lock() mpp = cfg.PARAMS['mp_processes'] - log.workflow('Initializing multiprocessing pool with ' - 'N={} processes.'.format(mpp)) _mp_pool = mp.Pool(mpp, initializer=_init_pool_globals, initargs=(cfg_contents, global_lock)) return _mp_pool @@ -219,7 +217,7 @@ def init_glacier_regions(rgidf=None, *, reset=False, force=False, prepro_rgi_version=None, prepro_base_url=None, from_tar=False, delete_tar=False, use_demo_glaciers=None): - """Initializes the list of Glacier Directories for this run. + """DEPRECATED: Initializes the list of Glacier Directories for this run. This is the very first task to do (always). If the directories are already available in the working directory, use them. If not, create new ones. @@ -264,6 +262,13 @@ def init_glacier_regions(rgidf=None, *, reset=False, force=False, ------- gdirs : list of :py:class:`oggm.GlacierDirectory` objects the initialised glacier directories + + Notes + ----- + This task is deprecated in favor of the more explicit + init_glacier_directories. Indeed, init_glacier_directories is very + similar to init_glacier_regions, but it does not process the DEMs: + a glacier directory is valid also without DEM. """ if reset and not force: @@ -283,9 +288,6 @@ def init_glacier_regions(rgidf=None, *, reset=False, force=False, if os.path.exists(fpath): rmtree(fpath) - # Read the hash dictionary before we use multiproc - utils.get_dl_verify_data('cluster.klima.uni-bremen.de') - gdirs = [] new_gdirs = [] if rgidf is None: @@ -296,8 +298,7 @@ def init_glacier_regions(rgidf=None, *, reset=False, force=False, # The dirs should be there already gl_dir = os.path.join(cfg.PATHS['working_dir'], 'per_glacier') for root, _, files in os.walk(gl_dir): - if files and ('outlines.shp' in files or - 'outlines.tar.gz' in files): + if files and ('dem.tif' in files): gdirs.append(oggm.GlacierDirectory(os.path.basename(root))) else: @@ -317,6 +318,8 @@ def init_glacier_regions(rgidf=None, *, reset=False, force=False, log.workflow('init_glacier_regions from prepro level {} on ' '{} glaciers.'.format(from_prepro_level, len(entities))) + # Read the hash dictionary before we use multiproc + utils.get_dl_verify_data('cluster.klima.uni-bremen.de') gdirs = execute_entity_task(gdir_from_prepro, entities, from_prepro_level=from_prepro_level, prepro_border=prepro_border, @@ -324,33 +327,168 @@ def init_glacier_regions(rgidf=None, *, reset=False, force=False, check_demo_glacier=use_demo_glaciers, base_url=prepro_base_url) else: - # TODO: if necessary this could use multiprocessing as well + # We can set the intersects file automatically here + if (cfg.PARAMS['use_intersects'] and + len(cfg.PARAMS['intersects_gdf']) == 0): + rgi_ids = np.unique(np.sort([entity.RGIId for entity in + entities])) + rgi_version = rgi_ids[0].split('-')[0][-2:] + fp = utils.get_rgi_intersects_entities(rgi_ids, + version=rgi_version) + cfg.set_intersects_db(fp) + for entity in entities: gdir = oggm.GlacierDirectory(entity, reset=reset, from_tar=from_tar, delete_tar=delete_tar) - outlines_path = gdir.get_filepath('outlines') - if not (os.path.exists(outlines_path) or - os.path.exists(outlines_path.replace('.shp', - '.tar.gz'))): - new_gdirs.append((gdir, dict(entity=entity))) + if not os.path.exists(gdir.get_filepath('dem')): + new_gdirs.append(gdir) gdirs.append(gdir) - # We can set the intersects file automatically here - if (cfg.PARAMS['use_intersects'] and new_gdirs and - (len(cfg.PARAMS['intersects_gdf']) == 0)): - rgi_ids = np.unique(np.sort([t[0].rgi_id for t in new_gdirs])) - rgi_version = new_gdirs[0][0].rgi_version - fp = utils.get_rgi_intersects_entities(rgi_ids, version=rgi_version) - cfg.set_intersects_db(fp) - if len(new_gdirs) > 0: + # Read the hash dictionary before we use multiproc + utils.get_dl_verify_data('cluster.klima.uni-bremen.de') # If not initialized, run the task in parallel execute_entity_task(tasks.define_glacier_region, new_gdirs) return gdirs +def init_glacier_directories(rgidf=None, *, reset=False, force=False, + from_prepro_level=None, prepro_border=None, + prepro_rgi_version=None, prepro_base_url=None, + from_tar=False, delete_tar=False, + use_demo_glaciers=None): + """Initializes the list of Glacier Directories for this run. + + This is the very first task to do (always). If the directories are already + available in the working directory, use them. If not, create new ones. + + Parameters + ---------- + rgidf : GeoDataFrame or list of ids, optional for pre-computed runs + the RGI glacier outlines. If unavailable, OGGM will parse the + information from the glacier directories found in the working + directory. It is required for new runs. + reset : bool + delete the existing glacier directories if found. + force : bool + setting `reset=True` will trigger a yes/no question to the user. Set + `force=True` to avoid this. + from_prepro_level : int + get the gdir data from the official pre-processed pool. See the + documentation for more information + prepro_border : int + for `from_prepro_level` only: if you want to override the default + behavior which is to use `cfg.PARAMS['border']` + prepro_rgi_version : str + for `from_prepro_level` only: if you want to override the default + behavior which is to use `cfg.PARAMS['rgi_version']` + prepro_base_url : str + for `from_prepro_level` only: if you want to override the default + URL from which to download the gdirs. Default currently is + https://cluster.klima.uni-bremen.de/~fmaussion/gdirs/oggm_v1.1/ + use_demo_glaciers : bool + whether to check the demo glaciers for download (faster than the + standard prepro downloads). The default is to decide whether or + not to check based on simple criteria such as glacier list size. + from_tar : bool, default=False + extract the gdir data from a tar file. If set to `True`, + will check for a tar file at the expected location in `base_dir`. + delete_tar : bool, default=False + delete the original tar file after extraction. + + Returns + ------- + gdirs : list of :py:class:`oggm.GlacierDirectory` objects + the initialised glacier directories + + Notes + ----- + This task is very similar to init_glacier_regions, with one main + difference: it does not process the DEMs for this glacier. + Eventually, init_glacier_regions will be deprecated and removed from the + codebase. + """ + + if reset and not force: + reset = utils.query_yes_no('Delete all glacier directories?') + + if prepro_border is None: + prepro_border = int(cfg.PARAMS['border']) + + if from_prepro_level and prepro_border not in [10, 80, 160, 250]: + if 'test' not in utils._downloads.GDIR_URL: + raise InvalidParamsError("prepro_border or cfg.PARAMS['border'] " + "should be one of: 10, 80, 160, 250.") + + # if reset delete also the log directory + if reset: + fpath = os.path.join(cfg.PATHS['working_dir'], 'log') + if os.path.exists(fpath): + rmtree(fpath) + + if rgidf is None: + # Infer the glacier directories from folders available in working dir + if reset: + raise ValueError('Cannot use reset without setting rgidf') + log.workflow('init_glacier_directories by parsing all available ' + 'folders (this takes time: if possible, provide rgidf ' + 'instead).') + # The dirs should be there already + gl_dir = os.path.join(cfg.PATHS['working_dir'], 'per_glacier') + gdirs = [] + for root, _, files in os.walk(gl_dir): + if files and ('outlines.shp' in files or + 'outlines.tar.gz' in files): + gdirs.append(oggm.GlacierDirectory(os.path.basename(root))) + else: + # Create glacier directories from input + # Check if dataframe or list of str + try: + entities = [] + for _, entity in rgidf.iterrows(): + entities.append(entity) + except AttributeError: + entities = utils.tolist(rgidf) + + # Check demo + if use_demo_glaciers is None: + use_demo_glaciers = len(entities) < 100 + + if from_prepro_level is not None: + log.workflow('init_glacier_directories from prepro level {} on ' + '{} glaciers.'.format(from_prepro_level, + len(entities))) + # Read the hash dictionary before we use multiproc + utils.get_dl_verify_data('cluster.klima.uni-bremen.de') + gdirs = execute_entity_task(gdir_from_prepro, entities, + from_prepro_level=from_prepro_level, + prepro_border=prepro_border, + prepro_rgi_version=prepro_rgi_version, + check_demo_glacier=use_demo_glaciers, + base_url=prepro_base_url) + else: + # We can set the intersects file automatically here + if (cfg.PARAMS['use_intersects'] and + len(cfg.PARAMS['intersects_gdf']) == 0): + rgi_ids = np.unique(np.sort([entity.RGIId for entity in + entities])) + rgi_version = rgi_ids[0].split('-')[0][-2:] + fp = utils.get_rgi_intersects_entities(rgi_ids, + version=rgi_version) + cfg.set_intersects_db(fp) + + gdirs = [] + for entity in entities: + gdir = oggm.GlacierDirectory(entity, reset=reset, + from_tar=from_tar, + delete_tar=delete_tar) + gdirs.append(gdir) + + return gdirs + + def gis_prepro_tasks(gdirs): """Shortcut function: run all flowline preprocessing tasks. @@ -361,6 +499,7 @@ def gis_prepro_tasks(gdirs): """ task_list = [ + tasks.define_glacier_region, tasks.glacier_masks, tasks.compute_centerlines, tasks.initialize_flowlines,