workflow/preprocessing/04_l2_select_and_zip.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select glaciers within region shapes and basic attributes: Level 2 files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L0 and L1 files are too large: they are selected roughly by lon lat boxes and a buffer, so there is a lot of duplication. Here we select glaciers by region by overlaying the glaciers representative points with the region outlines. This reduces the size of the files and make them ready for regional processing.\n",
    "\n",
    "In addition, we compute glacier area and store the CenLon, CenLat attributes for later use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import utils \n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import shutil\n",
    "import glob\n",
    "import os\n",
    "import subprocess\n",
    "import tarfile\n",
    "import shapely.geometry as shpg\n",
    "import progressbar\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def xy_coord(geom):\n",
    "    x, y = geom.xy\n",
    "    return x[0], y[0]"
   ]
  },
  {
   "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": [
    "reg_file = os.path.join(data_dir, 'l0_regions', '00_rgi70_regions', '00_rgi70_O1Regions.shp')\n",
    "reg_f = gpd.read_file(reg_file)\n",
    "# We drop the second alaska box which is useless\n",
    "reg_f = reg_f.drop(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idir = utils.mkdir(data_dir + '/l1_with_interiors')\n",
    "odir = utils.mkdir(data_dir + '/l2_sel_reg')\n",
    "odir_tar = utils.mkdir(data_dir + '/l2_sel_reg_tars')\n",
    "\n",
    "#for reg in range(1, 20):\n",
    "for reg in [10]:\n",
    "\n",
    "    fp = f'{idir}/RGI{reg:02d}/RGI{reg:02d}.shp'\n",
    "    print(fp)\n",
    "    print('Reading...')\n",
    "    shp = gpd.read_file(fp)\n",
    "    \n",
    "    print('Pointing...')\n",
    "    \n",
    "    rp = shp.representative_point()\n",
    "    \n",
    "    coordinates = np.array(list(rp.apply(xy_coord)))\n",
    "    shp['CenLon'] = coordinates[:, 0]\n",
    "    shp['CenLat'] = coordinates[:, 1]\n",
    "    \n",
    "    rp = rp.to_frame('geometry')\n",
    "    rp['orig_index'] = shp.index\n",
    "    \n",
    "    intersect = gpd.overlay(rp, reg_f.loc[reg_f.RGI_CODE == f'{reg:02d}'], how='intersection')\n",
    "    \n",
    "    odf = shp.loc[intersect['orig_index']]\n",
    "    \n",
    "    print('Area-ing...')\n",
    "    odf['area'] = odf.to_crs({'proj':'cea'}).area\n",
    "    \n",
    "    dd = utils.mkdir(f'{odir}/RGI{reg:02d}/', reset=True)\n",
    "    \n",
    "    print('Writing...')\n",
    "    odf.to_file(dd + f'RGI{reg:02d}.shp')\n",
    "    \n",
    "    print('Taring...')\n",
    "    print(subprocess.run(['tar', '-zcvf', f'{odir_tar}/RGI{reg:02d}.tar.gz', '-C', odir, f'RGI{reg:02d}']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}