A few more updates (#1726)

README.rst

@@ -69,18 +69,10 @@ About
         :target: https://github.com/OGGM/oggm/actions/workflows/run-tests.yml
         :alt: Linux build status
 
-    .. image:: https://img.shields.io/badge/Cross-validation-blue.svg
-        :target: https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/crossval.html
-        :alt: Mass balance cross validation
-
     .. image:: https://readthedocs.org/projects/oggm/badge/?version=latest
         :target: http://docs.oggm.org/en/latest/
         :alt: Documentation status
 
-    .. image:: https://img.shields.io/badge/benchmarked%20by-asv-green.svg?style=flat
-        :target: https://cluster.klima.uni-bremen.de/~github/asv/
-        :alt: Benchmark status
-
 :License:
     .. image:: https://img.shields.io/pypi/l/oggm.svg
         :target: https://github.com/OGGM/oggm/blob/master/LICENSE.txt

docs/environment.yml

@@ -2,7 +2,8 @@ name: oggm_docs
 channels:
   - conda-forge
 dependencies:
-  - numpy
+  - python
+  - numpy<2.0
   - scipy
   - pandas
   - shapely

docs/recommended_env.yml

@@ -2,7 +2,8 @@ name: oggm_env
 channels:
   - conda-forge
 dependencies:
-  - numpy
+  - python<3.12
+  - numpy<2.0
   - scipy
   - pandas
   - shapely
@@ -28,4 +29,5 @@ dependencies:
   - pip:
     - joblib
     - progressbar2
+    - git+https://github.com/OGGM/pytest-mpl
     - oggm

docs/recommended_minimal_env.yml

@@ -2,7 +2,7 @@ name: oggm_minimal
 channels:
   - conda-forge
 dependencies:
-  - numpy
+  - numpy<2.0
   - scipy
   - pandas
   - matplotlib

oggm/cli/benchmark.py

@@ -9,6 +9,7 @@
 import sys
 import argparse
 import time
+import platform
 import logging
 import pandas as pd
 import geopandas as gpd
@@ -73,7 +74,7 @@ def run_benchmark(rgi_version=None, rgi_reg=None, border=None,
     cfg.initialize(logging_level=logging_level, params=override_params)
 
     # Use multiprocessing?
-    cfg.PARAMS['use_multiprocessing'] = True
+    cfg.PARAMS['use_multiprocessing'] = platform.system() != 'Darwin'
 
     # How many grid points around the glacier?
     # Make it large if you expect your glaciers to grow large

oggm/core/centerlines.py

@@ -698,7 +698,7 @@ def _make_costgrid(mask, ext, z):
     # arbitrary high to avoid the lines to jump over adjacent boundaries
     cost[np.where(ext)] = np.nanmax(cost[np.where(ext)]) * 50
 
-    return np.where(mask, cost, np.Inf)
+    return np.where(mask, cost, np.inf)
 
 
 def _get_terminus_coord(gdir, ext_yx, zoutline):
@@ -1019,8 +1019,8 @@ def compute_downstream_line(gdir):
     _y = [ymesh[:, 0], ymesh[0, :], ymesh[:, -1], ymesh[-1, :]]
 
     # Find the way out of the domain
-    min_cost = np.Inf
-    min_len = np.Inf
+    min_cost = np.inf
+    min_len = np.inf
     line = None
     for h, x, y in zip(_h, _x, _y):
         ids = scipy.signal.argrelmin(h, order=10, mode='wrap')
@@ -1680,25 +1680,14 @@ def catchment_intersections(gdir):
 
     # We project them onto the mercator proj before writing. This is a bit
     # inefficient (they'll be projected back later), but it's more sustainable
-    try:
-        # salem for geopandas > 0.7
-        salem.transform_geopandas(gdfc, from_crs=gdir.grid,
-                                  to_crs=gdir.grid.proj, inplace=True)
-        salem.transform_geopandas(gdfi, from_crs=gdir.grid,
-                                  to_crs=gdir.grid.proj, inplace=True)
-    except TypeError:
-        # from_crs not available yet
-        if Version(gpd.__version__) >= Version('0.7.0'):
-            raise ImportError('You have installed geopandas v0.7 or higher. '
-                              'Please also update salem for compatibility.')
-        gdfc.crs = gdir.grid
-        gdfi.crs = gdir.grid
-        salem.transform_geopandas(gdfc, to_crs=gdir.grid.proj, inplace=True)
-        salem.transform_geopandas(gdfi, to_crs=gdir.grid.proj, inplace=True)
+    salem.transform_geopandas(gdfc, from_crs=gdir.grid,
+                              to_crs=gdir.grid.proj, inplace=True)
+    salem.transform_geopandas(gdfi, from_crs=gdir.grid,
+                              to_crs=gdir.grid.proj, inplace=True)
     if hasattr(gdfc.crs, 'srs'):
         # salem uses pyproj
-        gdfc.crs = gdfc.crs.srs
-        gdfi.crs = gdfi.crs.srs
+        gdfc.set_crs(crs=gdfc.crs.srs, inplace=True, allow_override=True)
+        gdfi.set_crs(crs=gdfi.crs.srs, inplace=True, allow_override=True)
     gdir.write_shapefile(gdfc, 'flowline_catchments')
     if len(gdfi) > 0:
         gdir.write_shapefile(gdfi, 'catchments_intersects')

oggm/core/gis.py

@@ -1994,9 +1994,11 @@ def reproject_gridded_data_variable_to_grid(gdir,
                            target_grid.dx ** 2)
 
             # only preserve total if there is some data before
-            factor = np.where(np.isclose(total_before, 0, atol=1e-6),
-                              0.,
-                              total_before / total_after)
+            with warnings.catch_warnings():
+                # Divide by zero is fine
+                warnings.filterwarnings("ignore", category=RuntimeWarning)
+                factor = np.where(np.isclose(total_before, 0, atol=1e-6),
+                                  0., total_before / total_after)
 
             if len(data.dims) == 3:
                 # need to add two axis for broadcasting

oggm/tests/test_prepro.py

@@ -755,6 +755,10 @@ def test_baltoro_centerlines(self):
         gdf['Status'] = '0'
         gdf['O1Region'] = '01'
         gdf['O2Region'] = '01'
+        del gdf['GLIMSID']
+        del gdf['CENLAT']
+        del gdf['CENLON']
+
         entity = gdf.iloc[0]
         gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir)
         gis.define_glacier_region(gdir)

oggm/tests/test_shop.py

@@ -10,7 +10,7 @@
 import pandas as pd
 from oggm import utils
 from oggm.utils import get_demo_file
-from oggm.shop import its_live, rgitopo, bedtopo, millan22, hugonnet_maps, glathida
+from oggm.shop import its_live, rgitopo, bedtopo, millan22, cook23, hugonnet_maps, glathida
 from oggm.core import gis, centerlines
 from oggm import cfg, tasks, workflow
 
@@ -166,6 +166,44 @@ def test_thickvel_to_glacier(self, class_case_dir, monkeypatch):
         assert df['millan_vol_km3'] > 174
 
 
+class Test_cook23:
+
+    @pytest.mark.slow
+    def test_cook23_to_glacier(self, class_case_dir, monkeypatch):
+
+        # Init
+        cfg.initialize()
+        cfg.PATHS['working_dir'] = class_case_dir
+        cfg.PARAMS['use_intersects'] = False
+        cfg.PATHS['dem_file'] = get_demo_file('hef_srtm.tif')
+        cfg.PARAMS['border'] = 10
+
+        entity = gpd.read_file(get_demo_file('Hintereisferner_RGI5.shp')).iloc[0]
+        gdir = oggm.GlacierDirectory(entity)
+        tasks.define_glacier_region(gdir)
+        tasks.glacier_masks(gdir)
+
+        # use our files
+        base_url = 'https://cluster.klima.uni-bremen.de/~oggm/test_files/cook23/'
+        monkeypatch.setattr(cook23, 'default_base_url', base_url)
+
+        cook23.cook23_to_gdir(gdir)
+
+        with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
+            mask = ds.glacier_mask.data.astype(bool)
+            thick = ds.cook23_thk.where(mask).data
+
+        # Simply some coverage and sanity checks
+        assert np.isfinite(thick).sum() / mask.sum() > 0.98
+        assert np.nanmax(thick) > 120
+        assert np.nansum(thick) * gdir.grid.dx**2 * 1e-9 > 0.4
+
+        # Stats
+        df = cook23.compile_cook23_statistics([gdir]).iloc[0]
+        assert df['cook23_perc_cov'] > 0.98
+        assert df['cook23_vol_km3'] > 0.4
+
+
 class Test_HugonnetMaps:
 
     @pytest.mark.slow
@@ -591,7 +629,7 @@ def test_all_at_once(self, class_case_dir):
         dfp = pd.concat(dfp, axis=1, keys=exps)
 
         # Common period
-        dfy = dft.resample('AS').mean().dropna().iloc[1:]
+        dfy = dft.resample('YS').mean().dropna().iloc[1:]
         dfm = dft.groupby(dft.index.month).mean()
         assert dfy.corr().min().min() > 0.44  # ERA5L and CERA do no correlate
         assert dfm.corr().min().min() > 0.97
@@ -610,7 +648,7 @@ def test_all_at_once(self, class_case_dir):
 
         # PRECIP
         # Common period
-        dfy = dfp.resample('AS').mean().dropna().iloc[1:] * 12
+        dfy = dfp.resample('YS').mean().dropna().iloc[1:] * 12
         dfm = dfp.groupby(dfp.index.month).mean()
         assert dfy.corr().min().min() > 0.5
         assert dfm.corr().min().min() > 0.8

oggm/tests/test_utils.py

@@ -4,6 +4,7 @@
 import time
 import hashlib
 import tarfile
+import platform
 import pytest
 import itertools
 from unittest import mock
@@ -876,6 +877,7 @@ def setUp(self):
         self.monkeypatch = MonkeyPatch()
         self.testdir = os.path.join(get_test_dir(), 'tmp_prepro_levs')
         self.reset_dir()
+        self.on_mac = platform.system() == 'Darwin'
 
     def tearDown(self):
         if os.path.exists(self.testdir):
@@ -1252,6 +1254,7 @@ def test_full_run_cru_centerlines(self):
                           mb_calibration_strategy='melt_temp',
                           test_intersects_file=inter,
                           test_topofile=topof,
+                          disable_mp=self.on_mac,
                           centerlines=True,
                           override_params={'geodetic_mb_period': ref_period,
                                            'baseline_climate': 'CRU',
@@ -1577,6 +1580,7 @@ def test_geodetic_per_glacier_and_massredis_run(self):
                           output_folder=odir, working_dir=wdir, is_test=True,
                           test_rgidf=rgidf, test_intersects_file=inter,
                           override_params=params,
+                          disable_mp=self.on_mac,
                           mb_calibration_strategy='melt_temp',
                           test_topofile=topof, elev_bands=True)
 
@@ -1671,6 +1675,7 @@ def test_start_from_prepro(self):
                           test_rgidf=rgidf, test_intersects_file=inter,
                           start_level=1, start_base_url=base_url,
                           mb_calibration_strategy='melt_temp',
+                          disable_mp=self.on_mac,
                           logging_level='INFO', max_level=5, elev_bands=True,
                           override_params={'geodetic_mb_period': ref_period,
                                            'baseline_climate': 'CRU',

oggm/tests/test_workflow.py

@@ -405,7 +405,7 @@ def test_merge_gridded_data():
     inv_volume_gridded_merged = (ds_merged.distributed_thickness.sum() *
                                  ds_merged.salem.grid.dx**2) * 1e-9
     assert_allclose(df['inv_volume_km3'].sum(), inv_volume_gridded_merged,
-                    rtol=2e-7)
+                    rtol=1e-6)
 
 
 @pytest.mark.slow

oggm/utils/_downloads.py

@@ -2018,7 +2018,7 @@ def get_rgi_glacier_entities(rgi_ids, version=None):
 
     # Make a new dataframe of those
     selection = pd.concat(selection)
-    selection.crs = sh.crs  # for geolocalisation
+    selection.set_crs(crs=sh.crs, inplace=True, allow_override=True)  # for geolocalisation
     if len(selection) != len(rgi_ids):
         raise RuntimeError('Could not find all RGI ids')
 
@@ -2226,7 +2226,7 @@ def get_rgi_intersects_entities(rgi_ids, version=None):
 
     # Make a new dataframe of those
     selection = pd.concat(selection)
-    selection.crs = sh.crs  # for geolocalisation
+    selection.set_crs(crs=sh.crs, inplace=True, allow_override=True)  # for geolocalisation
 
     return selection
 

oggm/utils/_funcs.py

@@ -532,8 +532,7 @@ def recursive_valid_polygons(geoms, crs):
         try:
             new_geoms.extend(recursive_valid_polygons(list(new_geom.geoms), crs))
         except AttributeError:
-            new_s = gpd.GeoSeries(new_geom)
-            new_s.crs = crs
+            new_s = gpd.GeoSeries(new_geom, crs=crs)
             if new_s.to_crs({'proj': 'cea'}).area.iloc[0] >= 10000:
                 new_geoms.append(new_geom)
     assert np.all([type(geom) == shpg.Polygon for geom in new_geoms])

oggm/utils/_workflow.py

@@ -784,7 +784,9 @@ def _write_shape_to_disk(gdf, fpath, to_tar=False):
     if '.shp' not in fpath:
         raise ValueError('File ending should be .shp')
 
-    gdf.to_file(fpath)
+    with warnings.catch_warnings():
+        warnings.filterwarnings('ignore', 'GeoSeries.notna', UserWarning)
+        gdf.to_file(fpath)
 
     if not to_tar:
         # Done here
@@ -890,9 +892,8 @@ def write_centerlines_to_shape(gdirs, *, path=True, to_tar=False,
                                 to_crs=_to_crs)
     # filter for none
     olist = [o for o in olist if o is not None]
-    odf = gpd.GeoDataFrame(itertools.chain.from_iterable(olist))
+    odf = gpd.GeoDataFrame(itertools.chain.from_iterable(olist), crs=to_crs)
     odf = odf.sort_values(by=['RGIID', 'SEGMENT_ID'])
-    odf.crs = to_crs
     # Sanity checks to avoid bad surprises
     gtype = np.array([g.geom_type for g in odf.geometry])
     if 'GeometryCollection' in gtype:
@@ -2461,7 +2462,7 @@ def idealized_gdir(surface_h, widths_m, map_dx, flowline_dx=1,
     entity.O1Region = '00'
     entity.O2Region = '0'
     gdir = GlacierDirectory(entity, base_dir=base_dir, reset=reset)
-    gdir.write_shapefile(gpd.GeoDataFrame([entity]), 'outlines')
+    gdir.write_shapefile(gpd.GeoDataFrame([entity], crs='EPSG:4326'), 'outlines')
 
     # Idealized flowline
     coords = np.arange(0, len(surface_h) - 0.5, 1)
@@ -2958,7 +2959,7 @@ def _reproject_and_write_shapefile(self, entity):
             if type(s) in [np.int32, np.int64]:
                 entity[k] = int(s)
         towrite = gpd.GeoDataFrame(entity).T.set_geometry('geometry')
-        towrite.crs = proj4_str
+        towrite.set_crs(crs=proj4_str, inplace=True, allow_override=True)
 
         # Write shapefile
         self.write_shapefile(towrite, 'outlines')