add region selection

bopen · Mar 19, 2024 · d25fd14 · d25fd14
1 parent bb1f7d5
commit d25fd14
Showing 1 changed file with 81 additions and 21 deletions.
diff --git a/notebooks/wp5/glacier_mass_balance_timeseries.ipynb b/notebooks/wp5/glacier_mass_balance_timeseries.ipynb
@@ -23,7 +23,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "import cartopy.crs as ccrs\n",
+    "import geopandas as gpd\n",
     "import matplotlib.pyplot as plt\n",
+    "import regionmask\n",
     "import scipy.stats\n",
     "import xarray as xr\n",
     "from c3s_eqc_automatic_quality_control import download\n",
@@ -48,15 +51,62 @@
    "source": [
     "period_start = \"1975_1976\"\n",
     "period_stop = \"2020_2021\"\n",
-    "assert all(\"_\" in period and len(period) == 9 for period in (period_start, period_stop))"
+    "assert all(\"_\" in period and len(period) == 9 for period in (period_start, period_stop))\n",
+    "\n",
+    "# Shapefile with regions\n",
+    "shapefile_url = \"https://www.gtn-g.ch/database/GlacReg_2017.zip\"\n",
+    "\n",
+    "# Select regions\n",
+    "regions_abbrev = [\n",
+    "    \"ALA\",\n",
+    "    \"CAU\",\n",
+    "    \"ASE\",\n",
+    "    \"ASC\",\n",
+    "    \"WNA\",\n",
+    "    \"SJM\",\n",
+    "    \"ACN\",\n",
+    "    \"ACS\",\n",
+    "    \"RUA\",\n",
+    "    \"ASN\",\n",
+    "    \"TRP\",\n",
+    "    \"GRL\",\n",
+    "    \"ISL\",\n",
+    "    \"CEU\",\n",
+    "    \"SAN\",\n",
+    "    \"ANT\",\n",
+    "    \"SCA\",\n",
+    "    \"ASW\",\n",
+    "    \"NZL\",\n",
+    "]\n",
+    "assert set(regions_abbrev) & {\n",
+    "    \"ALA\",\n",
+    "    \"CAU\",\n",
+    "    \"ASE\",\n",
+    "    \"ASC\",\n",
+    "    \"WNA\",\n",
+    "    \"SJM\",\n",
+    "    \"ACN\",\n",
+    "    \"ACS\",\n",
+    "    \"RUA\",\n",
+    "    \"ASN\",\n",
+    "    \"TRP\",\n",
+    "    \"GRL\",\n",
+    "    \"ISL\",\n",
+    "    \"CEU\",\n",
+    "    \"SAN\",\n",
+    "    \"ANT\",\n",
+    "    \"SCA\",\n",
+    "    \"ASW\",\n",
+    "    \"NZL\",\n",
+    "}"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "5",
    "metadata": {},
    "source": [
-    "## Define request"
+    "## Select regions"
    ]
   },
   {
@@ -66,26 +116,25 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "y0_start, y1_start = map(int, period_start.split(\"_\"))\n",
-    "y0_stop, y1_stop = map(int, period_stop.split(\"_\"))\n",
-    "collection_id = \"derived-gridded-glacier-mass-change\"\n",
-    "request = {\n",
-    "    \"variable\": \"glacier_mass_change\",\n",
-    "    \"product_version\": \"wgms_fog_2022_09\",\n",
-    "    \"format\": \"zip\",\n",
-    "    \"hydrological_year\": [\n",
-    "        f\"{y0}_{str(y1)[-2:]}\"\n",
-    "        for y0, y1 in zip(range(y0_start, y0_stop + 1), range(y1_start, y1_stop + 1))\n",
-    "    ],\n",
-    "}"
+    "gdf = gpd.read_file(shapefile_url)\n",
+    "gdf = gdf.dissolve(by=\"RGI_CODE\", as_index=False)\n",
+    "regions = regionmask.from_geopandas(\n",
+    "    gdf, numbers=\"RGI_CODE\", names=\"FULL_NAME\", abbrevs=\"WGMS_CODE\"\n",
+    ")[regions_abbrev]\n",
+    "ax = regions.plot(\n",
+    "    projection=ccrs.Robinson(),\n",
+    "    label=\"abbrev\",\n",
+    "    add_ocean=True,\n",
+    ")\n",
+    "ax.set_global()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "7",
    "metadata": {},
    "source": [
-    "## Function to cache"
+    "## Define request"
    ]
   },
   {
@@ -95,16 +144,26 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def spatial_sum(ds):\n",
-    "    return ds.sum((\"latitude\", \"longitude\"), keep_attrs=True)"
+    "y0_start, y1_start = map(int, period_start.split(\"_\"))\n",
+    "y0_stop, y1_stop = map(int, period_stop.split(\"_\"))\n",
+    "collection_id = \"derived-gridded-glacier-mass-change\"\n",
+    "request = {\n",
+    "    \"variable\": \"glacier_mass_change\",\n",
+    "    \"product_version\": \"wgms_fog_2022_09\",\n",
+    "    \"format\": \"zip\",\n",
+    "    \"hydrological_year\": [\n",
+    "        f\"{y0}_{str(y1)[-2:]}\"\n",
+    "        for y0, y1 in zip(range(y0_start, y0_stop + 1), range(y1_start, y1_stop + 1))\n",
+    "    ],\n",
+    "}"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "9",
    "metadata": {},
    "source": [
-    "## Download data"
+    "## Download and transform data"
    ]
   },
   {
@@ -118,8 +177,9 @@
     "    collection_id,\n",
     "    request,\n",
     "    chunks={\"hydrological_year\": 1},\n",
-    "    transform_func=spatial_sum,\n",
-    ").compute()\n",
+    ")\n",
+    "ds = ds.where(regions.mask(ds).notnull())\n",
+    "ds = ds.sum((\"latitude\", \"longitude\"), keep_attrs=True)\n",
     "\n",
     "# Customize attributes\n",
     "ds[\"time\"] = ds[\"time\"].dt.year\n",
@@ -271,7 +331,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.11.7"
   }
  },
  "nbformat": 4,