Change default map proj to utm (risky)

oggm/params.cfg

@@ -72,7 +72,7 @@ lru_maxsize = 100
 ### Glacier map and grid determination
 
 # Map projection: 'tmerc' or 'utm'
-map_proj = 'tmerc'
+map_proj = 'utm'
 
 # Decision on grid spatial resolution for each glacier
 # 'fixed': dx (meters) = fixed_dx

oggm/tests/funcs.py

@@ -409,6 +409,7 @@ def init_hef(reset=False, border=40, logging_level='INFO', rgi_id=None,
     cfg.PARAMS['prcp_fac'] = 2.5
     cfg.PARAMS['temp_bias_min'] = -10
     cfg.PARAMS['temp_bias_max'] = 10
+    cfg.PARAMS['map_proj'] = 'tmerc'
     hef_file = get_demo_file('Hintereisferner_RGI5.shp')
     entity = gpd.read_file(hef_file).iloc[0]
 
@@ -515,6 +516,7 @@ def init_columbia(reset=False):
     cfg.PARAMS['prcp_fac'] = 2.5
     cfg.PARAMS['baseline_climate'] = 'CRU'
     cfg.PARAMS['evolution_model'] = 'FluxBased'
+    cfg.PARAMS['map_proj'] = 'tmerc'
 
     entity = gpd.read_file(get_demo_file('01_rgi60_Columbia.shp')).iloc[0]
     gdir = oggm.GlacierDirectory(entity, reset=reset)
@@ -556,6 +558,7 @@ def init_columbia_eb(dir_name, reset=False):
     cfg.PARAMS['prcp_fac'] = 2.5
     cfg.PARAMS['baseline_climate'] = 'CRU'
     cfg.PARAMS['evolution_model'] = 'FluxBased'
+    cfg.PARAMS['map_proj'] = 'tmerc'
 
     entity = gpd.read_file(get_demo_file('01_rgi60_Columbia.shp')).iloc[0]
     gdir = oggm.GlacierDirectory(entity)

oggm/tests/test_benchmarks.py

@@ -383,7 +383,7 @@ def to_optimize(x):
             ds['ref'].data[df['j'], df['i']] = df['thick']
 
         rmsd = ((df.oggm - df.thick) ** 2).mean() ** .5
-        assert rmsd < 30
+        assert rmsd < 32
 
         dfm = df.mean()
         np.testing.assert_allclose(dfm.thick, dfm.oggm, 10)

oggm/tests/test_graphics.py

@@ -182,6 +182,7 @@ def test_multiple_inversion():
     cfg.PARAMS['use_winter_prcp_fac'] = False
     cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
+    cfg.PARAMS['map_proj'] = 'tmerc'
     cfg.PATHS['working_dir'] = testdir
 
     # Get the RGI ID
@@ -269,6 +270,7 @@ def test_multiple_models():
     cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
     cfg.PARAMS['border'] = 40
+    cfg.PARAMS['map_proj'] = 'tmerc'
 
     # Get the RGI ID
     hef_rgi = gpd.read_file(get_demo_file('divides_hef.shp'))
@@ -379,6 +381,7 @@ def test_chhota_shigri():
     cfg.PARAMS['use_winter_prcp_fac'] = False
     cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
+    cfg.PARAMS['map_proj'] = 'tmerc'
 
     hef_file = get_demo_file('divides_RGI50-14.15990.shp')
     df = gpd.read_file(hef_file)
@@ -424,6 +427,7 @@ def test_ice_cap():
     cfg.PARAMS['use_winter_prcp_fac'] = False
     cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
+    cfg.PARAMS['map_proj'] = 'tmerc'
 
     df = gpd.read_file(get_demo_file('divides_RGI50-05.08389.shp'))
     df['Area'] = df.Area * 1e-6  # cause it was in m2
@@ -467,6 +471,7 @@ def test_coxe():
     cfg.PARAMS['use_winter_prcp_fac'] = False
     cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
+    cfg.PARAMS['map_proj'] = 'tmerc'
 
     hef_file = get_demo_file('rgi_RGI50-01.10299.shp')
     entity = gpd.read_file(hef_file).iloc[0]

oggm/tests/test_prepro.py

@@ -69,6 +69,7 @@ def setUp(self):
         cfg.PATHS['dem_file'] = get_demo_file('hef_srtm.tif')
         cfg.PATHS['working_dir'] = self.testdir
         cfg.PARAMS['border'] = 20
+        cfg.PARAMS['use_multiprocessing'] = False
 
     def tearDown(self):
         self.rm_dir()
@@ -121,7 +122,7 @@ def test_define_region(self):
         # From string
         gdir = oggm.GlacierDirectory(gdir.rgi_id, base_dir=self.testdir)
         # This is not guaranteed to be equal because of projection issues
-        np.testing.assert_allclose(extent, gdir.extent_ll, atol=1e-5)
+        np.testing.assert_allclose(extent, gdir.extent_ll, atol=5e-4)
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", category=FutureWarning)
             # Warning in salem
@@ -263,7 +264,7 @@ def test_repr(self):
           Status: Glacier or ice cap
           Area: 8.036 km2
           Lon, Lat: (10.7584, 46.8003)
-          Grid (nx, ny): (159, 114)
+          Grid (nx, ny): (158, 116)
           Grid (dx, dy): (50.0, -50.0)
         """)
 
@@ -368,11 +369,12 @@ def test_compute_hypsometry_attributes(self):
 
         np.testing.assert_allclose(dfh['slope_deg'], entity.Slope, atol=0.5)
         np.testing.assert_allclose(dfh['aspect_deg'], entity.Aspect, atol=5)
-        np.testing.assert_allclose(dfh['zmed_m'], entity.Zmed, atol=20)
-        np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=20)
-        np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=20)
-        np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=20)
-        np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=20)
+        atol = 25
+        np.testing.assert_allclose(dfh['zmed_m'], entity.Zmed, atol=atol)
+        np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=atol)
+        np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=atol)
+        np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=atol)
+        np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=atol)
         # From google map checks
         np.testing.assert_allclose(dfh['terminus_lon'], 10.80, atol=0.01)
         np.testing.assert_allclose(dfh['terminus_lat'], 46.81, atol=0.01)
@@ -430,7 +432,7 @@ def test_rasterio_glacier_masks(self):
         gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir)
         gis.define_glacier_region(gdir)
 
-        # specifying a source will look for a DEN in a respective folder
+        # specifying a source will look for a DEM in a respective folder
         self.assertRaises(ValueError, gis.rasterio_glacier_mask,
                           gdir, source='SRTM')
 
@@ -461,7 +463,7 @@ def test_rasterio_glacier_masks(self):
 
         # not sure if we want such a hard coded test, but this will fail if the
         # sample data changes but could also indicate changes in rasterio
-        self.assertTrue(data.sum() == 3218)
+        np.testing.assert_allclose(data.sum(), 3218, atol=5)
 
     def test_intersects(self):
 
@@ -1166,7 +1168,7 @@ def test_width(self):
                          50.)
         h1, b = np.histogram(hgt, weights=harea, density=True, bins=bins)
         h2, b = np.histogram(rhgt, density=True, bins=bins)
-        assert utils.rmsd(h1*100*50, h2*100*50) < 1
+        assert utils.rmsd(h1*100*50, h2*100*50) < 1.1
 
         # Check that utility function is doing what is expected
         hh, ww = gdir.get_inversion_flowline_hw()
@@ -1921,8 +1923,8 @@ def to_optimize(x):
                                     tol=1e-4)['x']
         self.assertTrue(out[0] > 0.1)
         self.assertTrue(out[1] > 0.1)
-        self.assertTrue(out[0] < 1.1)
-        self.assertTrue(out[1] < 1.1)
+        self.assertTrue(out[0] < 1.2)
+        self.assertTrue(out[1] < 1.2)
         glen_a = glen_a * out[0]
         fs = fs * out[1]
         v = inversion.mass_conservation_inversion(gdir, fs=fs,
@@ -2004,7 +2006,7 @@ def test_invert_hef_from_consensus(self):
         df = workflow.calibrate_inversion_from_consensus(gdir,
                                                          a_bounds=a,
                                                          error_on_mismatch=False)
-        np.testing.assert_allclose(df.vol_itmix_m3, df.vol_oggm_m3, rtol=0.08)
+        np.testing.assert_allclose(df.vol_itmix_m3, df.vol_oggm_m3, rtol=0.1)
 
         # With fs it can work
         a = (0.1, 3)
@@ -2223,22 +2225,24 @@ def test_distribute(self):
         inversion.prepare_for_inversion(gdir)
 
         ref_v = 0.573 * 1e9
+        cfg.PARAMS['inversion_fs'] = 5.7e-20
 
         def to_optimize(x):
             glen_a = cfg.PARAMS['inversion_glen_a'] * x[0]
             fs = cfg.PARAMS['inversion_fs'] * x[1]
             v = inversion.mass_conservation_inversion(gdir, fs=fs,
-                                                         glen_a=glen_a)
+                                                      glen_a=glen_a)
             return (v - ref_v)**2
+
         import scipy.optimize as optimization
         out = optimization.minimize(to_optimize, [1, 1],
                                     bounds=((0.01, 10), (0.01, 10)),
                                     tol=1e-1)['x']
         glen_a = cfg.PARAMS['inversion_glen_a'] * out[0]
         fs = cfg.PARAMS['inversion_fs'] * out[1]
         v = inversion.mass_conservation_inversion(gdir, fs=fs,
-                                                     glen_a=glen_a,
-                                                     write=True)
+                                                  glen_a=glen_a,
+                                                  write=True)
         np.testing.assert_allclose(ref_v, v)
 
         inversion.distribute_thickness_interp(gdir, varname_suffix='_interp')
@@ -2287,7 +2291,7 @@ def to_optimize(x):
             glen_a = cfg.PARAMS['inversion_glen_a'] * x[0]
             fs = 0.
             v = inversion.mass_conservation_inversion(gdir, fs=fs,
-                                                         glen_a=glen_a)
+                                                      glen_a=glen_a)
             return (v - ref_v)**2
 
         import scipy.optimize as optimization
@@ -2296,13 +2300,13 @@ def to_optimize(x):
                                     tol=1e-4)['x']
 
         self.assertTrue(out[0] > 0.1)
-        self.assertTrue(out[0] < 10)
+        self.assertTrue(out[0] < 15)
 
         glen_a = cfg.PARAMS['inversion_glen_a'] * out[0]
         fs = 0.
         v = inversion.mass_conservation_inversion(gdir, fs=fs,
-                                                     glen_a=glen_a,
-                                                     write=True)
+                                                  glen_a=glen_a,
+                                                  write=True)
         np.testing.assert_allclose(ref_v, v)
 
         cls = gdir.read_pickle('inversion_output')

oggm/tests/test_shop.py

@@ -139,7 +139,7 @@ def test_thickvel_to_glacier(self, class_case_dir, monkeypatch):
             thick = ds.millan_ice_thickness.where(mask).data
 
         # Simply some coverage and sanity checks
-        assert np.isfinite(thick).sum() / mask.sum() > 0.92
+        assert np.isfinite(thick).sum() / mask.sum() > 0.97
         assert np.nanmax(thick) > 800
         assert np.nansum(thick) * gdir.grid.dx**2 * 1e-9 > 174
 
@@ -153,14 +153,14 @@ def test_thickvel_to_glacier(self, class_case_dir, monkeypatch):
             vy = ds.millan_vy.where(mask).data
 
         # Simply some coverage and sanity checks
-        assert np.isfinite(v).sum() / mask.sum() > 0.92
+        assert np.isfinite(v).sum() / mask.sum() > 0.97
         assert np.nanmax(v) > 2000
         assert np.nanmax(vx) > 500
         assert np.nanmax(vy) > 400
 
         # Stats
         df = millan22.compile_millan_statistics([gdir]).iloc[0]
-        assert df['millan_avg_vel'] > 180
+        assert df['millan_avg_vel'] > 175
         assert df['millan_max_vel'] > 2000
         assert df['millan_perc_cov'] > 0.95
         assert df['millan_vel_perc_cov'] > 0.92

oggm/tests/test_utils.py

@@ -1709,6 +1709,7 @@ def test_start_from_prepro(self):
                                            'baseline_climate': 'CRU',
                                            'use_winter_prcp_fac': False,
                                            'use_temp_bias_from_file': False,
+                                           'map_proj': 'tmerc',
                                            'prcp_fac': 2.5,
                                            }
                           )

oggm/tests/test_workflow.py

@@ -103,6 +103,7 @@ def up_to_climate(reset=False, use_mp=None, params_file=None):
     cfg.PARAMS['baseline_climate'] = 'CRU'
     cfg.PARAMS['evolution_model'] = 'FluxBased'
     cfg.PARAMS['downstream_line_shape'] = 'parabola'
+    cfg.PARAMS['map_proj'] = 'tmerc'
 
     # Go
     gdirs = workflow.init_glacier_directories(rgidf)

oggm/utils/_workflow.py

@@ -3424,7 +3424,18 @@ def _read_shapefile_from_path(cls, fp):
             if cfg.PARAMS['use_compression']:
                 fp += '.gz'
 
-        shp = gpd.read_file(fp)
+        retries = 5
+        delay = 0.1
+        for attempt in range(retries):
+            try:
+                shp = gpd.read_file(fp)
+                break
+            except (ValueError, IndexError, OSError) as err:
+                if attempt < retries - 1:
+                    logging.debug(f"Retrying read_file due to: {err} (attempt {attempt + 1})")
+                    time.sleep(delay)
+                else:
+                    raise
 
         # .properties file is created for compressed shapefiles. github: #904
         _properties = fp.replace('tar://', '') + '.properties'