add support of WGMS data in the geodetic scripts

Author: albangossard

massbalancemachine/data_preprocessing/wgms.py

@@ -0,0 +1,57 @@
+import os, sys
+import glob
+import pandas as pd
+import geopandas as gpd
+
+from data_processing.Dataset import Dataset
+from data_processing.utils import get_rgi
+from data_processing.wgms import wgms_folder
+from data_processing.glacier_utils import get_region_name
+
+processed_stakes_folder = os.path.join(wgms_folder, "processed")
+
+
+def processed_features_stakes_path(rgi_region):
+    if rgi_region is None:
+        return os.path.join(processed_stakes_folder, "all.csv")
+    else:
+        assert isinstance(rgi_region, int)
+        return os.path.join(processed_stakes_folder, f"region_{rgi_region}.csv")
+
+
+def build_monthly_data(data, cfg, rgi_region=None):
+
+    assert (
+        rgi_region is not None
+    ), "For the moment only one single RGI region can be used at a time with the WGMS data."
+
+    data = get_rgi(data=data, region=rgi_region)
+
+    # Drop measurements where no RGIId was found
+    data = data[data.RGIId.notnull()]
+
+    region_name = get_region_name(rgi_region)
+    dataset = Dataset(cfg, data=data, region_name=region_name, region_id=rgi_region)
+
+    voi_topographical = ["aspect", "slope", "svf"]
+
+    # Retrieve the topographical features for each stake measurement based on the latitude and longitude of the stake and add them to the dataset
+    dataset.get_topo_features(vois=voi_topographical)
+
+    df = dataset.data
+    df["MONTH_START"] = [str(date)[4:6] for date in df.FROM_DATE]
+    df["MONTH_END"] = [str(date)[4:6] for date in df.TO_DATE]
+    # df.MONTH_START.unique(), df.MONTH_END.unique()
+
+    dataset.get_climate_features()
+
+    # Specify the short names of the climate variables available in the dataset
+    vois_climate = ["t2m", "tp", "slhf", "sshf", "ssrd", "fal", "str"]
+
+    # For each record, convert to a monthly time resolution
+    dataset.convert_to_monthly(
+        vois_climate=vois_climate, vois_topographical=voi_topographical
+    )
+
+    os.makedirs(processed_stakes_folder, exist_ok=True)
+    dataset.data.to_csv(processed_features_stakes_path(rgi_region), index=False)

massbalancemachine/data_processing/glacier_utils.py

@@ -8,6 +8,7 @@
 import oggm
 import venv
 import subprocess
+from typing import Union
 
 import config
 from data_processing.product_utils import mbm_path
@@ -107,6 +108,15 @@ def get_region_shape_file(region: str):
     return shp_path
 
 
+def get_region_name(region: Union[str, int]):
+    if isinstance(region, int):
+        region = f"{region:02d}"
+    rgi_version = "62"
+    fp = oggm.utils.get_rgi_region_file(region, version=rgi_version)
+    rgi_name = os.path.basename(fp).split("_", 1)[1].replace(".shp", "").split("_")[1]
+    return rgi_name
+
+
 def get_region_area_bounds(region):
     if not isinstance(region, str):
         region = f"{region:02d}"

massbalancemachine/dataloader/SourceManager.py

@@ -26,6 +26,8 @@
 from data_processing.Product import Product
 from dataloader.DataLoader import DataLoader, set_dataloader_splits
 import data_preprocessing.iceland
+import data_preprocessing.wgms
+import data_processing.wgms
 import config
 
 ###
@@ -107,7 +109,12 @@
 )
 
 _default_test_glaciers_iceland = []
-_default_train_glaciers_iceland = ["RGI60-06.00228", "RGI60-06.00232"]
+_default_train_glaciers_iceland = [
+    "RGI60-06.00228",
+    "RGI60-06.00232",
+    "RGI60-06.00236",
+    "RGI60-06.00238",
+]
 
 _default_test_glaciers_norway = [
     "RGI60-08.01258",
@@ -178,6 +185,8 @@ def _default_input(sourceData):
         return []
     elif sourceData == "norway":
         return []
+    elif "wgms" in sourceData:
+        return []
     else:
         raise ValueError(f"source_data={sourceData} is unknown")
 
@@ -202,8 +211,8 @@ def __init__(self, cfg, params, test_split_on="GLACIER"):
         self.glaciers_list = []
         self.mean_stakes_elevation = {}
 
-    def train_test_glaciers(self):
-        return self.train_glaciers, self.test_glaciers
+    # def train_test_glaciers(self):
+    #     return self.train_glaciers, self.test_glaciers
 
     def load_stakes_data(self):
         raise NotImplemented(
@@ -236,12 +245,34 @@ def train_test_sets(self):
             meta_data_columns=self.cfg.metaData,
         )
 
-        # Ensure all test glaciers exist in the dataset
         existing_glaciers = set(
             dataloader.data.GLACIER.unique()
             if self.glacierNamesProvided
             else dataloader.data.RGIId.unique()
         )
+
+        if self.test_glaciers is None:
+            assert (
+                self.train_glaciers is not None
+            ), "If test_glaciers is not defined, train_glaciers should be defined to automatically determine the test glaciers based on available data."
+            self.test_glaciers = list(
+                set(existing_glaciers).difference(self.train_glaciers)
+            )
+            print(
+                f"Determining test glaciers based on available data: {self.test_glaciers}"
+            )
+        if self.train_glaciers is None:
+            assert (
+                self.test_glaciers is not None
+            ), "If train_glaciers is not defined, test_glaciers should be defined to automatically determine the train glaciers based on available data."
+            self.train_glaciers = list(
+                set(existing_glaciers).difference(self.test_glaciers)
+            )
+            print(
+                f"Determining train glaciers based on available data: {self.train_glaciers}"
+            )
+
+        # Ensure all test glaciers exist in the dataset
         missing_glaciers = [g for g in self.test_glaciers if g not in existing_glaciers]
 
         if missing_glaciers:
@@ -568,3 +599,47 @@ def load_stakes_data(self):
         data_monthly = dataloader.data
 
         return data_monthly
+
+
+class SourceManagerWGMS(SourceManager):
+    def __init__(self, cfg, params, *args, rgi_region=None, **kwargs):
+        self.train_glaciers = params["training"].get("train_glaciers")
+        self.test_glaciers = params["training"].get("test_glaciers")
+        self.rgi_region = rgi_region
+        super().__init__(cfg, params, *args, **kwargs)
+
+    def load_stakes_data(self):
+        # TODO: for the moment the arguments are the RGIId but we should manage this properly in the future
+
+        path_preprocessed = data_preprocessing.wgms.processed_features_stakes_path(
+            self.rgi_region
+        )
+        p = Product(path_preprocessed)
+        if not p.is_up_to_date():
+            data = data_processing.wgms.load_processed_wgms(rgi_region=self.rgi_region)
+            data_preprocessing.wgms.build_monthly_data(data, self.cfg, self.rgi_region)
+            p.gen_chk()
+        data = pd.read_csv(path_preprocessed)
+        _format_data_credit(
+            "WGMS (2026): Fluctuations of Glaciers (FoG) Database. World Glacier Monitoring Service (WGMS), Zurich, Switzerland. https://doi.org/10.5904/wgms-fog-2026-02-10",
+            "Open access under the requirement of correct citation",
+            [
+                "WGMS (2023): Global Glacier Change Bulletin No. 5 (2020-2021). Michael Zemp, Isabelle Gärtner-Roer, Samuel U. Nussbaumer, Ethan Z. Welty, Inés Dussaillant, and Jacqueline Bannwart (eds.), ISC (WDS) / IUGG (IACS) / UNEP / UNESCO / WMO, World Glacier Monitoring Service, Zurich, Switzerland, 134 pp. Based on database version https://doi.org/10.5904/wgms-fog-2023-09.",
+                "WGMS (2013): Glacier Mass Balance Bulletin No. 12 (2010-2011). Michael Zemp, Samuel U. Nussbaumer, Kathrin Naegeli, Isabelle Gärtner-Roer, Frank Paul, Martin Hoelzle, and Wilfried Haeberli (eds.), ICSU (WDS) / IUGG (IACS) / UNEP / UNESCO / WMO, World Glacier Monitoring Service, Zurich, Switzerland, 106 pp. Based on database version https://doi.org/10.5904/wgms-fog-2013-11.",
+                "WGMS (2012): Fluctuations of Glaciers 2005-2010 (Vol. X): Michael Zemp, Holger Frey, Isabelle Gärtner-Roer, Samuel U. Nussbaumer, Martin Hoelzle, Frank Paul, and Wilfried Haeberli (eds.), ICSU (WDS) / IUGG (IACS) / UNEP / UNESCO / WMO, World Glacier Monitoring Service, Zurich, Switzerland. Based on database version https://doi.org/10.5904/wgms-fog-2012-11.",
+            ],
+        )
+
+        data["aspect"] = 180 * data["aspect"] / np.pi
+        data["slope"] = 180 * data["slope"] / np.pi
+        data["t2m"] = data["t2m"] - 273.15
+
+        dataloader = DataLoader(
+            self.cfg,
+            data=data,
+            random_seed=self.cfg.seed,
+            meta_data_columns=self.cfg.metaData,
+        )
+        data_monthly = dataloader.data
+
+        return data_monthly

massbalancemachine/dataloader/__init__.py

@@ -4,5 +4,6 @@
     SourceManagerSwitzerland,
     SourceManagerIceland,
     SourceManagerNorway,
+    SourceManagerWGMS,
 )
 from dataloader.GeoDataLoader import GeoDataLoader

massbalancemachine/plots/perf_plots.py

@@ -502,8 +502,14 @@ def predVSTruthTimeSeries(
             scores_predVSTruth[k]["winter"] = scores_winter[k]
         if hasSummer:
             scores_predVSTruth[k]["summer"] = scores_summer[k]
+
+    # Drop entries for which PERIOD is not included in ["annual", "winter", "summer"]
+    grouped_ids_drop = grouped_ids[
+        grouped_ids.PERIOD.isin(["annual", "winter", "summer"])
+    ]
+
     predVSTruth(
-        grouped_ids,
+        grouped_ids_drop,
         ax1,
         scores=scores_predVSTruth,
         hue="PERIOD",

massbalancemachine/training/training.py

@@ -332,9 +332,15 @@ def assessOnTest(log_dir, model, geodataloader_test, light=False):
     mse_winter, rmse_winter, mae_winter, pearson_corr_winter, r2_winter, bias_winter = (
         scores(predWinter, targetWinter)
     )
-    mse_summer, rmse_summer, mae_summer, pearson_corr_summer, r2_summer, bias_summer = (
-        scores(predSummer, targetSummer)
-    )
+    if len(targetSummer) > 0:
+        (
+            mse_summer,
+            rmse_summer,
+            mae_summer,
+            pearson_corr_summer,
+            r2_summer,
+            bias_summer,
+        ) = scores(predSummer, targetSummer)
 
     # Make plots
     plot_pred_vs_obs(log_dir, targetAll, predAll, {"rmse": rmse, "mae": mae})
@@ -348,7 +354,7 @@ def assessOnTest(log_dir, model, geodataloader_test, light=False):
         plt.savefig(os.path.join(log_dir, "geodetic_test.png"))
         plt.close(fig)
 
-    return {
+    stats = {
         "mse": mse.item(),
         "rmse": rmse.item(),
         "mae": mae.item(),
@@ -367,13 +373,15 @@ def assessOnTest(log_dir, model, geodataloader_test, light=False):
         "pearson_winter": pearson_corr_winter.item(),
         "r2_winter": r2_winter.item(),
         "bias_winter": bias_winter.item(),
-        "mse_summer": mse_summer.item(),
-        "rmse_summer": rmse_summer.item(),
-        "mae_summer": mae_summer.item(),
-        "pearson_summer": pearson_corr_summer.item(),
-        "r2_summer": r2_summer.item(),
-        "bias_summer": bias_summer.item(),
     }
+    if len(targetSummer) > 0:
+        stats["mse_summer"] = mse_summer.item()
+        stats["rmse_summer"] = rmse_summer.item()
+        stats["mae_summer"] = mae_summer.item()
+        stats["pearson_summer"] = pearson_corr_summer.item()
+        stats["r2_summer"] = r2_summer.item()
+        stats["bias_summer"] = bias_summer.item()
+    return stats
 
 
 def plot_pred_vs_obs(log_dir, target, pred, scores):

scripts/common.py

@@ -47,6 +47,8 @@ def parseParams(params):
             "batch_size": batch_size,
             "weight_decay": weight_decay,
             "scalingStakes": scalingStakes,
+            "test_glaciers": params["training"].get("test_glaciers"),
+            "train_glaciers": params["training"].get("train_glaciers"),
         },
     }
 

scripts/geo/eval.py

@@ -103,6 +103,19 @@
         ],
         notMetaDataNotFeatures=["POINT_BALANCE", "svf"],
     )
+elif "wgms" in sourceData:
+    cfg = mbm.Config(
+        metaData=[
+            "RGIId",
+            "ID",
+            "N_MONTHS",
+            "MONTHS",
+            "PERIOD",
+            "YEAR",
+            "POINT_ELEVATION",
+        ],
+        notMetaDataNotFeatures=["POINT_BALANCE", "svf"],
+    )
 else:
     raise ValueError(f"source_data={sourceData} is unknown")
 
@@ -131,6 +144,15 @@
     datasetManager = mbm.dataloader.SourceManagerNorway(
         cfg, params, test_split_on=keyGlacier
     )
+elif "wgms" in sourceData:
+    _split = sourceData.split(":")
+    if len(_split) > 1:
+        rgi_region = int(_split[1])
+    else:
+        rgi_region = None
+    datasetManager = mbm.dataloader.SourceManagerWGMS(
+        cfg, params, test_split_on="RGIId", rgi_region=rgi_region
+    )
 train_set, test_set, months_head_pad, months_tail_pad = datasetManager.train_test_sets()
 
 
@@ -173,6 +195,8 @@
         test_glaciers = list(data_test.GLACIER.unique())
     elif sourceData in ["iceland", "norway"]:
         test_glaciers = list(data_test.RGIId.unique())
+    elif "wgms" in sourceData:
+        test_glaciers = list(data_test.RGIId.unique())
 
     assert set(df_X_test_subset.RGIId.unique()) == set(test_glaciers)
 
@@ -204,11 +228,14 @@
         "r2": scores["winter"]["r2"],
         "bias": scores["winter"]["bias"],
     }
-    scores_summer = {
-        "rmse": scores["summer"]["rmse"],
-        "r2": scores["summer"]["r2"],
-        "bias": scores["summer"]["bias"],
-    }
+    if "summer" in scores:
+        scores_summer = {
+            "rmse": scores["summer"]["rmse"],
+            "r2": scores["summer"]["r2"],
+            "bias": scores["summer"]["bias"],
+        }
+    else:
+        scores_summer = None
 
     fig = mbm.plots.predVSTruthTimeSeries(
         grouped_ids=grouped_ids,
@@ -289,6 +316,8 @@
     train_glaciers = list(data_train.GLACIER.unique())
 elif sourceData in ["iceland", "norway"]:
     train_glaciers = list(data_train.RGIId.unique())
+elif "wgms" in sourceData:
+    train_glaciers = list(data_train.RGIId.unique())
 
 # Create dataloader
 train_gdl = mbm.dataloader.GeoDataLoader(
@@ -315,11 +344,14 @@
     "r2": scores_train["winter"]["r2"],
     "bias": scores_train["winter"]["bias"],
 }
-scores_summer = {
-    "rmse": scores_train["summer"]["rmse"],
-    "r2": scores_train["summer"]["r2"],
-    "bias": scores_train["summer"]["bias"],
-}
+if "summer" in scores_train:
+    scores_summer = {
+        "rmse": scores_train["summer"]["rmse"],
+        "r2": scores_train["summer"]["r2"],
+        "bias": scores_train["summer"]["bias"],
+    }
+else:
+    scores_summer = None
 
 fig = mbm.plots.predVSTruthTimeSeries(
     grouped_ids=grouped_ids_train,