workflow/preprocessing/01_rgi7_reg_files.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modify the RGI6 regions files for RGI7 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**List of changes from RGI6 to RGI7**: \n",
    "- The southern boundary of region 12 (Caucasus and Middle East) has been shifted south by 2° (from 32°N to 30°N) to encompass a cluster of glaciers which were previously unassigned.\n",
    "- The region boxes for region 01 (Alaska) used to encompass some islands in the Bering Sea East of Kamtchatka. This subregion contains no glaciers (and probably hasn't for a long time) and is now removed. \n",
    "- The data type of the `rgi_code` attribute in the first-order region file is now `str` (was `int`). The `rgi_code` now has a leading zero, for example `02` instead of `2`. \n",
    "- Region 19 has been renamed to \"*Subantarctic and Antarctic Islands*\" (was: \"*Antarctic and Subantarctic*\")\n",
    "- The polygon of RGI region 19 has been updated to remove antarctic mainland. Potentially, a new region (`20`, \"*Antarctic Mainland*\") might be used in the future if glaciers are inventorized there.\n",
    "- All abbreviations in the second-order regions file have been replaced by their full name (e.g. \"East Central\" instead of \"EC\")\n",
    "- The first-order and second-order region files now have a field called `long_code` which contains a string representing the full region name, using the lowercase with underscores format (e.g. `02_western_canada_usa`). This field is used to name the corresponding RGI shapefiles.\n",
    "- The subregion `05-11` (Greenland Ice Sheet) has been removed.\n",
    "- The `WGMS_CODE` column has been deleted from all files.\n",
    "- The `RGI_CODE` column is now called `o1region` (first order files) and `o2region` (second-order files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# go down from rgi7_scripts/workflow/preprocessing\n",
    "data_dir = '../../../rgi7_data/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import shapely.geometry as shpg\n",
    "import geopandas as gpd\n",
    "from utils import mkdir\n",
    "import shutil\n",
    "from shapely.validation import make_valid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_dir = os.path.abspath(os.path.join(data_dir, '00_rgi70_regions'))\n",
    "mkdir(out_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the new names\n",
    "df_ref = pd.read_csv(os.path.join(out_dir, '00_rgi70_rgi_region_names.csv'), dtype=str)\n",
    "df_ref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the RGI region files\n",
    "rgi_dir = os.path.join(data_dir, 'l0_RGIv6')\n",
    "rgi_reg = gpd.read_file('zip://' + os.path.join(data_dir, 'l0_RGIv6', '00_rgi60_regions.zip', '00_rgi60_O1Regions.shp'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new region 20 based on the subregion outline\n",
    "rgi_subreg = gpd.read_file('zip://' + os.path.join(data_dir, 'l0_RGIv6', '00_rgi60_regions.zip', '00_rgi60_O2Regions.shp'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_subreg.loc[rgi_subreg['RGI_CODE'] == '19-31'].plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aa_geom = rgi_subreg.loc[rgi_subreg['RGI_CODE'] == '19-31'].iloc[0].geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suba_geom = rgi_reg.iloc[-1].geometry.difference(aa_geom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(suba_geom.geoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suba_geom.geoms[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suba_geom = suba_geom.geoms[0]\n",
    "suba_geom.is_valid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suba_geom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg.loc[20, 'geometry'] = suba_geom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aa_geom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg = pd.concat([rgi_reg, gpd.GeoDataFrame([aa_geom], columns=['geometry'], crs=rgi_reg.crs)], ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(~ rgi_reg.is_valid).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg.loc[~ rgi_reg.is_valid]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, s in rgi_reg.loc[~ rgi_reg.is_valid].iterrows():\n",
    "    corr = make_valid(s.geometry)\n",
    "    if type(corr) == shpg.Polygon:\n",
    "        print(f'Success : {s.FULL_NAME}')\n",
    "        rgi_reg.loc[i, 'geometry'] = corr\n",
    "    else:\n",
    "        for g in corr.geoms:\n",
    "            print(type(g))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg.loc[~ rgi_reg.is_valid]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = make_valid(rgi_reg.loc[~ rgi_reg.is_valid].iloc[0].geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr.geoms[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr.geoms[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg.loc[21, 'geometry'] = corr.geoms[0]\n",
    "rgi_reg.loc[21, 'RGI_CODE'] = 20\n",
    "rgi_reg['RGI_CODE'] = rgi_reg['RGI_CODE'].astype(int)\n",
    "rgi_reg.loc[21, 'FULL_NAME'] = 'Antarctic Mainland'\n",
    "rgi_reg.loc[20, 'FULL_NAME'] = 'Subantarctic and Antarctic Islands'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove the useless Alaska box in the Eastern hemisphere\n",
    "rgi_reg = rgi_reg.drop(1).reset_index()\n",
    "del rgi_reg['index']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename things that need renaming\n",
    "rgi_reg['RGI_CODE'] = [f'{int(i):02d}' for i in rgi_reg['RGI_CODE']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg = rgi_reg.rename({'WGMS_CODE':'LONG_CODE'}, axis=1)\n",
    "rgi_reg['FULL_NAME'] = [df_ref.loc[df_ref.RGI_CODE == i, 'FULL_NAME'].values[0] for i in rgi_reg['RGI_CODE']]\n",
    "rgi_reg['LONG_CODE'] = [df_ref.loc[df_ref.RGI_CODE == i, 'LONG_CODE'].values[0] for i in rgi_reg['RGI_CODE']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the RGI 12 polygon\n",
    "poly = rgi_reg.loc[rgi_reg.RGI_CODE == '12'].iloc[0].geometry\n",
    "poly.bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go down to 30° South instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = poly.exterior.xy\n",
    "ny = np.where(np.isclose(y, 31), 30, y)\n",
    "new_poly = shpg.Polygon(np.array((x, ny)).T)\n",
    "rgi_reg.loc[rgi_reg.RGI_CODE == '12', 'geometry'] = new_poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(~ rgi_reg.is_valid).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove RGI 20\n",
    "# rgi_reg = rgi_reg.iloc[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remame columns\n",
    "names = [s.lower() for s in rgi_reg.columns]\n",
    "names[1] = 'o1region'\n",
    "rgi_reg.columns = names\n",
    "rgi_reg = rgi_reg[['o1region', 'full_name', 'long_code', 'geometry']]\n",
    "rgi_reg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(rgi_reg.set_index('o1region')[['full_name', 'long_code']].to_markdown())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg.crs = 'EPSG:4326'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_file_dir = os.path.join(out_dir, '00_rgi70_O1Regions')\n",
    "rgi_reg.to_file(out_file_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shutil.make_archive(f'{out_dir}/00_rgi70_O1Regions', 'zip', out_file_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check\n",
    "rgi_reg = gpd.read_file(out_file_dir)\n",
    "assert rgi_reg.o1region.dtype == 'O'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gpd.io.file.infer_schema(rgi_reg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subregions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg = gpd.read_file('zip://' + os.path.join(data_dir, 'l0_RGIv6', '00_rgi60_regions.zip', '00_rgi60_O2Regions.shp'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the new names\n",
    "df_ref = pd.read_csv(os.path.join(out_dir, '00_rgi70_rgi_subregion_names.csv'), dtype=str)\n",
    "df_ref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(rgi_reg.RGI_CODE.unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a long code\n",
    "long = []\n",
    "for s in df_ref['FULL_NAME']:\n",
    "    s = s.replace('(', '').replace(')', '').replace('/', '_').replace('-', '_')\n",
    "    s = '_'.join([so.lower() for so in s.split(' ')])\n",
    "    long.append(s.replace('__', '_').replace('__', '_'))\n",
    "df_ref['LONG_CODE'] = long"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Move RGI12 a bit south\n",
    "poly = rgi_reg.loc[rgi_reg.RGI_CODE == '12-02'].iloc[0].geometry\n",
    "poly.bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = poly.exterior.xy\n",
    "ny = np.where(np.isclose(y, 32), 30, y)\n",
    "new_poly = shpg.Polygon(np.array((x, ny)).T)\n",
    "rgi_reg.loc[rgi_reg.RGI_CODE == '12-02', 'geometry'] = new_poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove the useless Alaska box in the Eastern hemisphere\n",
    "rgi_reg = rgi_reg.drop(3).reset_index()\n",
    "del rgi_reg['index']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename things that need renaming\n",
    "rgi_reg.loc[91, 'RGI_CODE'] = '20-01'\n",
    "rgi_reg.loc[91, 'FULL_NAME'] = 'Antarctic Mainland'\n",
    "rgi_reg = rgi_reg.rename({'WGMS_CODE':'LONG_CODE'}, axis=1)\n",
    "rgi_reg['FULL_NAME'] = [df_ref.loc[df_ref.RGI_CODE == i, 'FULL_NAME'].values[0] for i in rgi_reg['RGI_CODE']]\n",
    "rgi_reg['LONG_CODE'] = [i + '_' + df_ref.loc[df_ref.RGI_CODE == i, 'LONG_CODE'].values[0] for i in rgi_reg['RGI_CODE']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(~ rgi_reg.is_valid).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, s in rgi_reg.loc[~ rgi_reg.is_valid].iterrows():\n",
    "    corr = make_valid(s.geometry)\n",
    "    if type(corr) == shpg.Polygon:\n",
    "        print(f'Success : {s.FULL_NAME}')\n",
    "        rgi_reg.loc[i, 'geometry'] = corr\n",
    "    else:\n",
    "        print(f'Not good : {s.FULL_NAME}')\n",
    "        for g in corr.geoms:\n",
    "            print(type(g))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr.geoms[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr.geoms[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg.loc[i, 'geometry'] = corr.geoms[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(~ rgi_reg.is_valid).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove RGI 20\n",
    "# rgi_reg = rgi_reg.iloc[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove Greenland icesheet\n",
    "rgi_reg = rgi_reg.loc[rgi_reg.RGI_CODE != '05-11'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add o1region for help\n",
    "rgi_reg['o1Region'] = [s.split('-')[0] for s in rgi_reg.RGI_CODE]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# order \n",
    "rgi_reg = rgi_reg[['o1Region', 'RGI_CODE', 'FULL_NAME', 'LONG_CODE', 'geometry']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remame things\n",
    "names = [s.lower() for s in rgi_reg.columns]\n",
    "names[1] = 'o2region'\n",
    "rgi_reg.columns = names\n",
    "rgi_reg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(rgi_reg.o2region.unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg.crs = 'EPSG:4326'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_file_dir = os.path.join(out_dir, '00_rgi70_O2Regions')\n",
    "rgi_reg.to_file(out_file_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shutil.make_archive(f'{out_dir}/00_rgi70_O2Regions', 'zip', out_file_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(rgi_reg.set_index('o2region')[['full_name', 'long_code']].to_markdown())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}