add bedmachine to prepro

OGGM · fmaussion · Nov 20, 2024 · Nov 20, 2024 · Nov 21, 2024 · Nov 21, 2024
commit d2b53d77a83a565a2091e124d697217e94d44cc5
diff --git a/docs/api.rst b/docs/api.rst
@@ -100,6 +100,8 @@ OGGM Shop
     shop.millan22.thickness_to_gdir
     shop.millan22.velocity_to_gdir
     shop.millan22.compile_millan_statistics
+    shop.bedmachine.bedmachine_to_gdir
+    shop.bedmachine.compile_bedmachine_statistics
     shop.rgitopo.init_glacier_directories_from_rgitopo
     shop.rgitopo.select_dem_from_dir
     shop.rgitopo.dem_quality_check

diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py
@@ -84,7 +84,8 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
                       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_glathida=False,
+                      add_hugonnet_dhdt=False, add_bedmachine=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,
@@ -159,6 +160,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_bedmachine : bool
+        adds (reprojects) the bedmachine ice thickness 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.
@@ -481,6 +485,10 @@ 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_bedmachine:
+            from oggm.shop.bedmachine import bedmachine_to_gdir
+            workflow.execute_entity_task(bedmachine_to_gdir, gdirs)
+            bin_variables.append('bedmachine_ice_thickness')
         if add_glathida:
             from oggm.shop.glathida import glathida_to_gdir
             workflow.execute_entity_task(glathida_to_gdir, gdirs)
@@ -550,6 +558,10 @@ def _time_log():
             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_bedmachine:
+            from oggm.shop.bedmachine import compile_bedmachine_statistics
+            opath = os.path.join(sum_dir, 'bedmachine_statistics_{}.csv'.format(rgi_reg))
+            compile_bedmachine_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))
@@ -899,7 +911,12 @@ def parse_args(args):
                              'With --elev-bands, the data will also be '
                              'binned.')
     parser.add_argument('--add-hugonnet-dhdt', nargs='?', const=True, default=False,
-                        help='adds (reprojects) the millan dhdt '
+                        help='adds (reprojects) the hugonnet dhdt '
+                             'maps to the glacier directories. '
+                             'With --elev-bands, the data will also be '
+                             'binned.')
+    parser.add_argument('--add-bedmachine', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the Bedmachine ice thickness '
                              'maps to the glacier directories. '
                              'With --elev-bands, the data will also be '
                              'binned.')
@@ -995,6 +1012,7 @@ def parse_args(args):
                 add_itslive_velocity=args.add_itslive_velocity,
                 add_millan_velocity=args.add_millan_velocity,
                 add_hugonnet_dhdt=args.add_hugonnet_dhdt,
+                add_bedmachine=args.add_bedmachine,
                 add_glathida=args.add_glathida,
                 dynamic_spinup=dynamic_spinup,
                 err_dmdtda_scaling_factor=args.err_dmdtda_scaling_factor,

diff --git a/oggm/shop/bedmachine.py b/oggm/shop/bedmachine.py
@@ -68,7 +68,7 @@ def bedmachine_to_gdir(gdir):
     # Write
     with utils.ncDataset(gdir.get_filepath('gridded_data'), 'a') as nc:
 
-        vn = 'bedmachine_thickness'
+        vn = 'bedmachine_ice_thickness'
         if vn in nc.variables:
             v = nc.variables[vn]
         else:
@@ -98,7 +98,7 @@ def bedmachine_statistics(gdir):
 
     try:
         with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
-            thick = ds['bedmachine_thickness'].where(ds['glacier_mask'], np.nan).load()
+            thick = ds['bedmachine_ice_thickness'].where(ds['glacier_mask'], np.nan).load()
             gridded_area = ds['glacier_mask'].sum() * gdir.grid.dx ** 2 * 1e-6
             d['bedmachine_area_km2'] = float((~thick.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6)
             d['bedmachine_perc_cov'] = float(d['bedmachine_area_km2'] / gridded_area)

diff --git a/oggm/shop/its_live.py b/oggm/shop/its_live.py
@@ -27,7 +27,7 @@
 regions = [f'RGI{reg:02d}A' for reg in range(1, 20)]
 
 rgi_region_links = {'01': 'RGI01A',
-                    '02': ['RGI01A', 'RGI02A'],
+                    '02': 'RGI02A',
                     '03': 'RGI03A',
                     '04': 'RGI04A',
                     '05': 'RGI05A',
@@ -46,7 +46,12 @@
 
 
 def _find_region(gdir):
-    return rgi_region_links.get(gdir.rgi_region, None)
+    reg_n = gdir.rgi_region
+    if reg_n == '02':
+        # Northermost glaciers are in reg 1 file
+        if gdir.cenlat > 55:
+            reg_n = '01'
+    return rgi_region_links.get(reg_n, None)
 
 
 def _region_files():