atlxi_lake.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "# **ICESat-2 Active Subglacial Lakes in Antarctica**\n",
    "\n",
    "Finding subglacial lakes that are draining or filling under the ice!\n",
    "They can be detected with ICESat-2 data, as significant changes in height\n",
    "(> 1 metre) over a relatively short duration (< 1 year), i.e. a high rate of\n",
    "elevation change over time (dhdt).\n",
    "\n",
    "In this notebook, we'll use some neat tools to help us examine the lakes:\n",
    "- To find active subglacial lake boundaries,\n",
    "use an *unsupervised clustering* technique\n",
    "- To see ice surface elevation trends at a higher temporal resolution (< 3 months),\n",
    "perform *crossover track error analysis* on intersecting ICESat-2 tracks\n",
    "\n",
    "To speed up analysis on millions of points,\n",
    "we will use state of the art GPU algorithms enabled by RAPIDS AI libraries,\n",
    "or parallelize the processing across our HPC's many CPU cores using Dask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import subprocess\n",
    "\n",
    "# os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"1\"\n",
    "\n",
    "import cudf\n",
    "import cuml\n",
    "import dask\n",
    "import dask.array\n",
    "import geopandas as gpd\n",
    "import hvplot.xarray\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import panel as pn\n",
    "import pygmt\n",
    "import scipy.spatial\n",
    "import shapely.geometry\n",
    "import tqdm\n",
    "import xarray as xr\n",
    "import zarr\n",
    "\n",
    "import deepicedrain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "use_cupy: bool = True\n",
    "if use_cupy:\n",
    "    import cupy\n",
    "    import dask_cuda\n",
    "    import dask_cudf\n",
    "\n",
    "    cluster = dask_cuda.LocalCUDACluster(n_workers=2, threads_per_worker=1)\n",
    "else:\n",
    "    cluster = dask.distributed.LocalCluster(n_workers=8, threads_per_worker=1)\n",
    "\n",
    "client = dask.distributed.Client(address=cluster)\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(\"ATLXI/df_dhdt_antarctica.parquet\"):\n",
    "    zarrarray = zarr.open_consolidated(store=f\"ATLXI/ds_dhdt_antarctica.zarr\", mode=\"r\")\n",
    "    _ = deepicedrain.ndarray_to_parquet(\n",
    "        ndarray=zarrarray,\n",
    "        parquetpath=\"ATLXI/df_dhdt_antarctica.parquet\",\n",
    "        variables=[\"x\", \"y\", \"dhdt_slope\", \"referencegroundtrack\"],  # \"h_corr\"],\n",
    "        dropnacols=[\"dhdt_slope\"],\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Read in Antarctic Drainage Basin Boundaries shapefile into a GeoDataFrame\n",
    "ice_boundaries: gpd.GeoDataFrame = (\n",
    "    deepicedrain.catalog.measures_antarctic_boundaries.read()\n",
    ")\n",
    "drainage_basins: gpd.GeoDataFrame = ice_boundaries.query(expr=\"TYPE == 'GR'\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load in ICESat-2 data (x, y, dhdt) and do initial trimming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trimmed 290925594 -> 90113189\n"
     ]
    }
   ],
   "source": [
    "# Read in raw x, y, dhdt_slope and referencegroundtrack data into the GPU\n",
    "cudf_raw: cudf.DataFrame = dask_cudf.read_parquet(\n",
    "    path=\"ATLXI/df_dhdt_antarctica.parquet\",\n",
    "    columns=[\"x\", \"y\", \"dhdt_slope\", \"referencegroundtrack\"],\n",
    "    # filters=[[('dhdt_slope', '<', -0.105)], [('dhdt_slope', '>', 0.105)]],\n",
    ")\n",
    "# Filter to points with dhdt that is less than -0.105 m/yr or more than +0.105 m/yr\n",
    "# Based on ICESat-2 ATL06's accuracy and precision of 3.3 ± 7.2cm from Brunt et al 2020\n",
    "# See https://doi.org/10.1029/2020GL090572\n",
    "cudf_many = cudf_raw.loc[abs(cudf_raw.dhdt_slope) > 0.105].compute()\n",
    "print(f\"Trimmed {len(cudf_raw)} -> {len(cudf_many)}\")\n",
    "if \"cudf_raw\" in globals():\n",
    "    del cudf_raw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clip outlier values to 3 sigma (standard deviations) from mean\n",
    "_mean = cudf_many.dhdt_slope.mean()\n",
    "_std = cudf_many.dhdt_slope.std()\n",
    "cudf_many.dhdt_slope.clip(\n",
    "    lower=np.float32(_mean - 3 * _std), upper=np.float32(_mean + 3 * _std), inplace=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Label ICESat-2 points according to their drainage basin\n",
    "\n",
    "Uses Point in Polygon.\n",
    "For each point, find out which Antarctic Drainage Basin they are in.\n",
    "This will also remove the points on floating (FR) ice shelves and islands (IS),\n",
    "so that we keep only points on the grounded (GR) ice regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trimmed 90113189 -> 55596512\n"
     ]
    }
   ],
   "source": [
    "# Use point in polygon to label points according to the drainage basins they fall in\n",
    "cudf_many[\"drainage_basin\"]: cudf.Series = deepicedrain.point_in_polygon_gpu(\n",
    "    points_df=cudf_many, poly_df=drainage_basins, poly_limit=16\n",
    ")\n",
    "X_many = cudf_many.dropna()  # drop points that are not in a drainage basin\n",
    "print(f\"Trimmed {len(cudf_many)} -> {len(X_many)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "# Find Active Subglacial Lake clusters\n",
    "\n",
    "Uses Density-based spatial clustering of applications with noise (DBSCAN)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subglacial Lake Finder algorithm\n",
    "\n",
    "For each Antarctic drainage basin:\n",
    "\n",
    "1. Select all points with significant elevation change over time (dhdt)\n",
    "  - Specifically, the (absolute) dhdt value should be\n",
    "    2x the median (absolute) dhdt for that drainage basin\n",
    "  - E.g. if median dhdt for basin is 0.35 m/yr,\n",
    "    we choose points that have dhdt > 0.70 m/yr\n",
    "2. Run unsupervised clustering to pick out active subglacial lakes\n",
    "  - Split into draining (-dhdt) and filling (+dhdt) points first\n",
    "  - Use DBSCAN algorithm to cluster points into groups,\n",
    "    with an eps (distance) of 3 km and minimum sample size of 250 points\n",
    "3. Check each potential point cluster to see if it meets active lake criteria\n",
    "  1. Build a convex hull 'lake' polygon around clustered points\n",
    "  2. Check that the 'lake' has significant elevation change relative to outside\n",
    "    - For the area in the 5 km buffer region **outside** the 'lake' polygon:\n",
    "       - Find median dhdt (outer_dhdt)\n",
    "       - Find median absolute deviation of dhdt values (outer_mad)\n",
    "    - For the area **inside** the 'lake' polygon:\n",
    "       - Find median dhdt (inner_dhdt)\n",
    "    - If the potential lake shows an elevation change that is more than\n",
    "      3x the surrounding deviation of background elevation change,\n",
    "      we infer that this is likely an active subglacial 'lake'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/198 [00:01<02:57,  1.10it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "102075 rows at Academy above ± 0.44 m/yr\n",
      "2 draining and 9 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  6%|▌         | 12/198 [00:02<00:25,  7.38it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "37918 rows at Jutulstraumen above ± 0.58 m/yr\n",
      "0 draining and 1 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 27%|██▋       | 54/198 [00:07<00:17,  8.03it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "70286 rows at Cook above ± 0.51 m/yr\n",
      "0 draining and 1 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 28%|██▊       | 56/198 [00:08<00:18,  7.64it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39403 rows at David above ± 0.50 m/yr\n",
      "1 draining and 1 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 30%|███       | 60/198 [00:09<00:31,  4.37it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90734 rows at Mercer above ± 0.55 m/yr\n",
      "5 draining and 15 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 31%|███▏      | 62/198 [00:09<00:24,  5.58it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "288050 rows at Pine_Island above ± 0.97 m/yr\n",
      "2 draining and 1 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 33%|███▎      | 66/198 [00:19<02:36,  1.19s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "160978 rows at Thwaites above ± 0.78 m/yr\n",
      "4 draining and 3 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 35%|███▌      | 70/198 [00:20<01:19,  1.62it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "126226 rows at Whillans above ± 0.64 m/yr\n",
      "6 draining and 13 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 36%|███▋      | 72/198 [00:23<01:50,  1.14it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "63649 rows at Kamb above ± 0.54 m/yr\n",
      "2 draining and 12 filling lakes found\n",
      "6238 rows at Leverett above ± 0.61 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 37%|███▋      | 74/198 [00:24<01:36,  1.28it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "86214 rows at Scott above ± 0.51 m/yr\n",
      "5 draining and 8 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 39%|███▉      | 77/198 [00:25<01:02,  1.93it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "58941 rows at Amundsen above ± 0.51 m/yr\n",
      "4 draining and 5 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 41%|████      | 81/198 [00:25<00:38,  3.01it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "77517 rows at Beardmore above ± 0.48 m/yr\n",
      "2 draining and 2 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 41%|████▏     | 82/198 [00:26<00:40,  2.84it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "66146 rows at Nimrod above ± 0.48 m/yr\n",
      "2 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 43%|████▎     | 86/198 [00:27<00:31,  3.61it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "96005 rows at Byrd above ± 0.47 m/yr\n",
      "5 draining and 5 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 45%|████▍     | 89/198 [00:27<00:23,  4.62it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26337 rows at Bindschadler above ± 0.46 m/yr\n",
      "2 draining and 1 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 46%|████▌     | 91/198 [00:28<00:24,  4.44it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "49745 rows at MacAyeal above ± 0.47 m/yr\n",
      "4 draining and 3 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 56%|█████▌    | 110/198 [00:29<00:07, 12.29it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7570 rows at Bailey above ± 0.49 m/yr\n",
      "0 draining and 1 filling lakes found\n",
      "47112 rows at Slessor above ± 0.47 m/yr\n",
      "12 draining and 5 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 57%|█████▋    | 112/198 [00:30<00:18,  4.60it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30240 rows at Support_Force above ± 0.43 m/yr\n",
      "3 draining and 2 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 58%|█████▊    | 115/198 [00:31<00:22,  3.62it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "96788 rows at Foundation above ± 0.43 m/yr\n",
      "2 draining and 12 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 59%|█████▉    | 117/198 [00:32<00:18,  4.45it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20534 rows at Lambert above ± 0.41 m/yr\n",
      "0 draining and 1 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 61%|██████    | 120/198 [00:32<00:13,  5.57it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19521 rows at Mellor above ± 0.48 m/yr\n",
      "1 draining and 1 filling lakes found\n",
      "9735 rows at Fisher above ± 0.49 m/yr\n",
      "0 draining and 1 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 67%|██████▋   | 133/198 [00:33<00:06, 10.45it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9630 rows at Moller above ± 0.46 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 68%|██████▊   | 135/198 [00:34<00:10,  6.14it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "31449 rows at Institute above ± 0.47 m/yr\n",
      "7 draining and 3 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 72%|███████▏  | 142/198 [00:35<00:09,  5.61it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30329 rows at Bowman_Strom_Live_Axel-Heigerg above ± 0.77 m/yr\n",
      "2 draining and 0 filling lakes found\n",
      "27320 rows at Sulzberger above ± 0.82 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 74%|███████▎  | 146/198 [00:36<00:16,  3.20it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "77210 rows at Getz above ± 1.55 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 82%|████████▏ | 162/198 [00:39<00:07,  4.55it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64469 rows at Recovery above ± 0.43 m/yr\n",
      "4 draining and 6 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 198/198 [00:42<00:00,  4.69it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total of 193 subglacial lakes found\n"
     ]
    }
   ],
   "source": [
    "# Subglacial lake finder\n",
    "activelakes: dict = {\n",
    "    \"basin_name\": [],  # Antarctic drainage basin name\n",
    "    \"refgtracks\": [],  # Pipe-delimited list of ICESat-2 reference ground tracks\n",
    "    \"num_points\": [],  # Number of clustered data points\n",
    "    \"maxabsdhdt\": [],  # Maximum absolute dhdt value inside of lake boundary\n",
    "    \"inner_dhdt\": [],  # Median elev change over time (dhdt) inside of lake bounds\n",
    "    \"mean_dhdt\": [],  # Mean elev change over time (dhdt) inside of lake bounds\n",
    "    \"outer_dhdt\": [],  # Median elevation change over time (dhdt) outside of lake\n",
    "    \"outer_std\": [],  # Standard deviation of dhdt outside of lake\n",
    "    \"outer_mad\": [],  # Median absolute deviation of dhdt outside of lake\n",
    "    \"geometry\": [],  # Shapely Polygon geometry holding lake boundary coordinates\n",
    "}\n",
    "basin_name: str = \"Cook\"  # Set a basin name here\n",
    "basins = drainage_basins[drainage_basins.NAME == basin_name].index  # one specific basin\n",
    "# basins = drainage_basins[\n",
    "#     drainage_basins.NAME.isin((\"Cook\", \"Whillans\"))\n",
    "# ].index  # some specific basins\n",
    "basins: pd.core.indexes.numeric.Int64Index = drainage_basins.index  # run on all basins\n",
    "\n",
    "eps: int = 3000  # ICESat-2 tracks are separated by ~3 km across track, with each laser pair ~90 m apart\n",
    "min_samples: int = 300\n",
    "for basin_index in tqdm.tqdm(iterable=basins):\n",
    "    # Initial data cleaning, filter to rows that are in the drainage basin\n",
    "    basin = drainage_basins.loc[basin_index]\n",
    "    X_local = X_many.loc[X_many.drainage_basin == basin.NAME]  # .reset_index(drop=True)\n",
    "\n",
    "    # Get points with dhdt_slope higher than 3x the median dhdt_slope for the basin\n",
    "    # E.g. if median dhdt_slope is 0.30 m/yr, then we cluster points over 0.90 m/yr\n",
    "    abs_dhdt = X_local.dhdt_slope.abs()\n",
    "    tolerance: float = 3 * abs_dhdt.median()\n",
    "    X = X_local.loc[abs_dhdt > tolerance]\n",
    "\n",
    "    if len(X) <= 1000:  # don't run on too few points\n",
    "        continue\n",
    "\n",
    "    # Run unsupervised clustering separately on draining and filling lakes\n",
    "    # Draining lake points have negative labels (e.g. -1, -2, 3),\n",
    "    # Filling lake points have positive labels (e.g. 1, 2, 3),\n",
    "    # Noise points have NaN labels (i.e. NaN)\n",
    "    cluster_vars = [\"x\", \"y\", \"dhdt_slope\"]\n",
    "    draining_lake_labels = -deepicedrain.find_clusters(\n",
    "        X=X.loc[X.dhdt_slope < 0][cluster_vars],\n",
    "        eps=eps,\n",
    "        min_samples=min_samples,\n",
    "        verbose=cuml.common.logger.level_error,\n",
    "    )\n",
    "    filling_lake_labels = deepicedrain.find_clusters(\n",
    "        X=X.loc[X.dhdt_slope > 0][cluster_vars],\n",
    "        eps=eps,\n",
    "        min_samples=min_samples,\n",
    "        verbose=cuml.common.logger.level_error,\n",
    "    )\n",
    "    lake_labels = cudf.concat(objs=[draining_lake_labels, filling_lake_labels])\n",
    "    lake_labels: cudf.Series = lake_labels.sort_index()\n",
    "    assert lake_labels.name == \"cluster_id\"\n",
    "\n",
    "    # Checking all potential subglacial lakes in a basin\n",
    "    clusters: cudf.Series = lake_labels.unique()\n",
    "    for cluster_label in clusters.to_array():\n",
    "        # Store attribute and geometry information of each active lake\n",
    "        lake_points: cudf.DataFrame = X.loc[lake_labels == cluster_label]\n",
    "\n",
    "        # More data cleaning, dropping clusters with too few points\n",
    "        try:\n",
    "            assert len(lake_points) > 100\n",
    "        except AssertionError:\n",
    "            lake_labels = lake_labels.replace(to_replace=cluster_label, value=None)\n",
    "            continue\n",
    "\n",
    "        multipoint: shapely.geometry.MultiPoint = shapely.geometry.MultiPoint(\n",
    "            points=lake_points[[\"x\", \"y\"]].as_matrix()\n",
    "        )\n",
    "        convexhull: shapely.geometry.Polygon = multipoint.convex_hull\n",
    "\n",
    "        # Filter out (most) false positive subglacial lakes\n",
    "        # Check that elevation change over time in lake is anomalous to outside\n",
    "        # The 5000 m distance from lake boundary setting is empirically based on\n",
    "        # Smith et al. 2009's methodology at https://doi.org/10.3189/002214309789470879\n",
    "        outer_ring_buffer = convexhull.buffer(distance=5000) - convexhull\n",
    "        X_local[\"in_donut_ring\"] = deepicedrain.point_in_polygon_gpu(\n",
    "            points_df=X_local,\n",
    "            poly_df=gpd.GeoDataFrame({\"name\": True, \"geometry\": [outer_ring_buffer]}),\n",
    "        )\n",
    "        outer_points = X_local.dropna(subset=\"in_donut_ring\")\n",
    "        outer_dhdt: float = outer_points.dhdt_slope.median()\n",
    "\n",
    "        outer_std: float = outer_points.dhdt_slope.std()\n",
    "        outer_mad: float = scipy.stats.median_abs_deviation(\n",
    "            x=outer_points.dhdt_slope.to_pandas()\n",
    "        )\n",
    "\n",
    "        mean_dhdt: float = lake_points.dhdt_slope.mean()\n",
    "        inner_dhdt: float = lake_points.dhdt_slope.median()\n",
    "        X_local = X_local.drop(labels=\"in_donut_ring\", axis=\"columns\")\n",
    "\n",
    "        # If lake interior's median dhdt value is within 3 median absolute deviations\n",
    "        # of the lake exterior's dhdt value, we remove the lake label\n",
    "        # I.e. skip if above background change not significant enough\n",
    "        # Inspired by Kim et al. 2016's methodology at https://doi.org/10.5194/tc-10-2971-2016\n",
    "        if abs(inner_dhdt - outer_dhdt) < 3 * outer_mad:\n",
    "            lake_labels = lake_labels.replace(to_replace=cluster_label, value=None)\n",
    "            continue\n",
    "\n",
    "        maxabsdhdt: float = (\n",
    "            lake_points.dhdt_slope.max()\n",
    "            if cluster_label > 0  # positive label = filling\n",
    "            else lake_points.dhdt_slope.min()  # negative label = draining\n",
    "        )\n",
    "        refgtracks: str = \"|\".join(\n",
    "            map(str, lake_points.referencegroundtrack.unique().to_pandas())\n",
    "        )\n",
    "\n",
    "        # Save key variables to dictionary that will later go into geodataframe\n",
    "        activelakes[\"basin_name\"].append(basin.NAME)\n",
    "        activelakes[\"refgtracks\"].append(refgtracks)\n",
    "        activelakes[\"num_points\"].append(len(lake_points))\n",
    "        activelakes[\"maxabsdhdt\"].append(maxabsdhdt)\n",
    "        activelakes[\"inner_dhdt\"].append(inner_dhdt)\n",
    "        activelakes[\"mean_dhdt\"].append(mean_dhdt)\n",
    "        activelakes[\"outer_dhdt\"].append(outer_dhdt)\n",
    "        activelakes[\"outer_std\"].append(outer_std)\n",
    "        activelakes[\"outer_mad\"].append(outer_mad)\n",
    "        activelakes[\"geometry\"].append(convexhull)\n",
    "\n",
    "    # Calculate total number of lakes found for one drainage basin\n",
    "    clusters: cudf.Series = lake_labels.unique()\n",
    "    n_draining, n_filling = (clusters < 0).sum(), (clusters > 0).sum()\n",
    "    if n_draining + n_filling > 0:\n",
    "        print(f\"{len(X)} rows at {basin.NAME} above ± {tolerance:.2f} m/yr\")\n",
    "        print(f\"{n_draining} draining and {n_filling} filling lakes found\")\n",
    "\n",
    "if len(activelakes[\"geometry\"]) >= 1:\n",
    "    gdf = gpd.GeoDataFrame(activelakes, crs=\"EPSG:3031\")\n",
    "    basename = \"antarctic_subglacial_lakes\"  # f\"temp_{basin_name.lower()}_lakes\"  #\n",
    "    gdf.to_file(filename=f\"{basename}_3031.geojson\", driver=\"GeoJSON\")\n",
    "    gdf.to_crs(crs={\"init\": \"epsg:4326\"}).to_file(\n",
    "        filename=f\"{basename}_4326.geojson\", driver=\"GeoJSON\"\n",
    "    )\n",
    "\n",
    "print(f\"Total of {len(gdf)} subglacial lakes found\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize lakes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Concatenate XY points with labels, and move data from GPU to CPU\n",
    "X: cudf.DataFrame = cudf.concat(objs=[X, lake_labels], axis=\"columns\")\n",
    "X_ = X.to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Plot clusters on a map in colour, noise points/outliers as small dots\n",
    "fig = pygmt.Figure()\n",
    "n_clusters_ = len(X_.cluster_id.unique()) - 1  # No. of clusters minus noise (NaN)\n",
    "sizes = (X_.cluster_id.isna()).map(arg={True: 0.01, False: 0.1})\n",
    "pygmt.makecpt(cmap=\"polar\", series=(-1, 1, 2), color_model=\"+cDrain,Fill\", reverse=True)\n",
    "fig.plot(\n",
    "    x=X_.x,\n",
    "    y=X_.y,\n",
    "    sizes=sizes,\n",
    "    style=\"cc\",\n",
    "    color=pd.cut(x=X_.cluster_id, bins=(-np.inf, 0, np.inf), labels=[-1, 1]),\n",
    "    cmap=True,\n",
    "    frame=[\n",
    "        f'WSne+t\"Estimated number of lake clusters at {basin.NAME}: {n_clusters_}\"',\n",
    "        'xafg+l\"Polar Stereographic X (m)\"',\n",
    "        'yafg+l\"Polar Stereographic Y (m)\"',\n",
    "    ],\n",
    ")\n",
    "basinx, basiny = basin.geometry.exterior.coords.xy\n",
    "fig.plot(x=basinx, y=basiny, pen=\"thinnest,-\")\n",
    "fig.colorbar(position='JMR+w2c/0.5c+m+n\"Unclassified\"', L=\"i0.5c\")\n",
    "fig.savefig(fname=f\"figures/subglacial_lake_clusters_at_{basin.NAME}.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Select a subglacial lake to examine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Load dhdt data from Parquet file\n",
    "placename: str = \"siple_coast\"  # \"slessor_downstream\"  #  \"Recovery\"  # \"Whillans\"\n",
    "df_dhdt: cudf.DataFrame = cudf.read_parquet(\n",
    "    f\"ATLXI/df_dhdt_{placename.lower()}.parquet\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "index                                                         0\n",
      "geometry      POLYGON ((-444731.6953220846 -545129.683759524...\n",
      "basin_name                                             Whillans\n",
      "refgtracks    74|135|196|266|327|388|577|638|769|830|1019|10...\n",
      "num_points                                                 3422\n",
      "maxabsdhdt                                             6.731061\n",
      "inner_dhdt                                             1.152791\n",
      "mean_dhdt                                              1.365484\n",
      "outer_dhdt                                             0.338404\n",
      "outer_std                                              0.151085\n",
      "outer_mad                                              0.081393\n",
      "Name: 0, dtype: object\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"300\" height=\"300\" viewBox=\"-456630.37936665915 -545831.6143030439 18952.12467501691 17100.473795304657\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,-1074562.754810783)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"126.34749783344606\" opacity=\"0.6\" d=\"M -444731.6953220846,-545129.6837595248 L -445011.25885498157,-545106.2968520958 L -445376.93950812094,-545031.6591196847 L -447409.3526933879,-544382.2805507262 L -450352.5856256729,-542856.9347974786 L -451613.65405471186,-542088.3567467475 L -453168.7111508558,-541005.8354607922 L -454541.14624465484,-540030.3540272998 L -455129.4508626142,-539161.6296595022 L -455928.44882314,-537936.6574710397 L -455889.53593075665,-537286.6599730677 L -454243.5090539303,-535519.7292632827 L -453242.79123223276,-534521.465268024 L -449104.8354354054,-531478.2867125623 L -448620.22966874443,-531145.8302439799 L -443146.6610664141,-529433.0710512584 L -441727.20288979204,-529844.5068711725 L -439855.5902259358,-530889.0817186871 L -439736.6727377982,-531015.2193078512 L -439630.14850606647,-531141.0021021602 L -438656.6297759172,-534025.5166154493 L -438380.1852351614,-535865.5798414922 L -439494.1268296778,-541257.5824335428 L -441276.22981477395,-542729.7436845967 L -444731.6953220846,-545129.6837595248 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x7fbc201cd670>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Choose one Antarctic active subglacial lake polygon with EPSG:3031 coordinates\n",
    "lake_name: str = \"Whillans IX\"\n",
    "lake_catalog = deepicedrain.catalog.subglacial_lakes()\n",
    "lake_ids, transect_id = (\n",
    "    pd.json_normalize(lake_catalog.metadata[\"lakedict\"])\n",
    "    .query(\"lakename == @lake_name\")[[\"ids\", \"transect\"]]\n",
    "    .iloc[0]\n",
    ")\n",
    "lake = (\n",
    "    lake_catalog.read()\n",
    "    .loc[lake_ids]\n",
    "    .dissolve(by=np.zeros(shape=len(lake_ids), dtype=\"int64\"), as_index=False)\n",
    "    .squeeze()\n",
    ")\n",
    "\n",
    "region = deepicedrain.Region.from_gdf(gdf=lake, name=lake_name)\n",
    "draining: bool = lake.inner_dhdt < 0\n",
    "\n",
    "print(lake)\n",
    "lake.geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 21/21 [00:00<00:00, 36.83it/s]\n"
     ]
    }
   ],
   "source": [
    "# Subset data to lake of interest\n",
    "placename: str = region.name.lower().replace(\" \", \"_\")\n",
    "df_lake: cudf.DataFrame = region.subset(data=df_dhdt)\n",
    "# Get all raw xyz points and one transect line dataframe\n",
    "track_dict: dict = deepicedrain.split_tracks(df=df_lake.to_pandas())\n",
    "track_points: pd.DataFrame = (\n",
    "    pd.concat(track_dict.values())\n",
    "    .groupby(by=[\"x\", \"y\"])\n",
    "    .mean()  # z value is mean h_corr over all cycles\n",
    "    .reset_index()[[\"x\", \"y\", \"h_corr\"]]\n",
    ")\n",
    "try:\n",
    "    _rgt, _pt = transect_id.split(\"_\")\n",
    "    df_transect: pd.DataFrame = (\n",
    "        track_dict[transect_id][[\"x\", \"y\", \"h_corr\", \"cycle_number\"]]\n",
    "        .groupby(by=[\"x\", \"y\"])\n",
    "        .max()  # z value is maximum h_corr over all cycles\n",
    "        .reset_index()\n",
    "    )\n",
    "except AttributeError:\n",
    "    pass\n",
    "\n",
    "# Save lake outline to OGR GMT file format\n",
    "outline_points: str = f\"figures/{placename}/{placename}.gmt\"\n",
    "if not os.path.exists(path=outline_points):\n",
    "    os.makedirs(name=f\"figures/{placename}\", exist_ok=True)\n",
    "    lake_catalog.read().loc[list(lake_ids)].to_file(\n",
    "        filename=outline_points, driver=\"OGR_GMT\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an interpolated ice surface elevation grid for each ICESat-2 cycle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 6/6 [00:03<00:00,  1.71it/s]\n"
     ]
    }
   ],
   "source": [
    "# Generate gridded time-series of ice elevation over lake\n",
    "cycles: tuple = (3, 4, 5, 6, 7, 8, 9)\n",
    "os.makedirs(name=f\"figures/{placename}\", exist_ok=True)\n",
    "ds_lake: xr.Dataset = deepicedrain.spatiotemporal_cube(\n",
    "    table=df_lake.to_pandas(),\n",
    "    placename=placename,\n",
    "    cycles=cycles,\n",
    "    folder=f\"figures/{placename}\",\n",
    ")\n",
    "ds_lake.to_netcdf(path=f\"figures/{placename}/xyht_{placename}.nc\", mode=\"w\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elevation limits are: (273.2418518066406, 310.9927062988281)\n"
     ]
    }
   ],
   "source": [
    "# Get 3D grid_region (xmin/xmax/ymin/ymax/zmin/zmax),\n",
    "# and calculate normalized z-values as Elevation delta relative to Cycle 3\n",
    "z_limits: tuple = (float(ds_lake.z.min()), float(ds_lake.z.max()))  # original z limits\n",
    "grid_region: tuple = region.bounds() + z_limits\n",
    "\n",
    "ds_lake_diff: xr.Dataset = ds_lake - ds_lake.sel(cycle_number=3).z\n",
    "z_diff_limits: tuple = (float(ds_lake_diff.z.min()), float(ds_lake_diff.z.max()))\n",
    "diff_grid_region: np.ndarray = np.append(arr=grid_region[:4], values=z_diff_limits)\n",
    "\n",
    "print(f\"Elevation limits are: {z_limits}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3D plot of mean ice surface elevation (<z>) and rate of ice elevation change (dhdt)\n",
    "fig = deepicedrain.plot_icesurface(\n",
    "    grid=f\"figures/{placename}/xyht_{placename}.nc?z_mean\",\n",
    "    grid_region=grid_region,\n",
    "    diff_grid=f\"figures/{placename}/xyht_{placename}.nc?dhdt\",\n",
    "    track_points=track_points.to_numpy(),\n",
    "    outline_points=outline_points,\n",
    "    azimuth=157.5,  # 202.5  # 270\n",
    "    elevation=45,  # 60\n",
    "    title=f\"{region.name} Ice Surface\",\n",
    ")\n",
    "# Plot crossing transect line\n",
    "fig.plot3d(\n",
    "    data=df_transect[[\"x\", \"y\", \"h_corr\"]].to_numpy(),\n",
    "    color=\"yellow2\",\n",
    "    style=\"c0.1c\",\n",
    "    zscale=True,\n",
    "    perspective=True,\n",
    ")\n",
    "fig.savefig(f\"figures/{placename}/dsm_{placename}_cycles_{cycles[0]}-{cycles[-1]}.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# 3D plots of gridded ice surface elevation over time (one per cycle)\n",
    "for cycle in tqdm.tqdm(iterable=cycles):\n",
    "    time_nsec: pd.Timestamp = df_lake[f\"utc_time_{cycle}\"].to_pandas().mean()\n",
    "    time_sec: str = np.datetime_as_string(arr=time_nsec.to_datetime64(), unit=\"s\")\n",
    "\n",
    "    # grid = ds_lake.sel(cycle_number=cycle).z\n",
    "    fig = deepicedrain.plot_icesurface(\n",
    "        grid=f\"figures/{placename}/h_corr_{placename}_cycle_{cycle}.nc\",\n",
    "        grid_region=grid_region,\n",
    "        diff_grid=ds_lake_diff.sel(cycle_number=cycle).z,\n",
    "        diff_grid_region=diff_grid_region,\n",
    "        track_points=df_lake[[\"x\", \"y\", f\"h_corr_{cycle}\"]].dropna().as_matrix(),\n",
    "        outline_points=outline_points,\n",
    "        azimuth=157.5,  # 202.5  # 270\n",
    "        elevation=45,  # 60\n",
    "        title=f\"{region.name} at Cycle {cycle} ({time_sec})\",\n",
    "    )\n",
    "    fig.savefig(f\"figures/{placename}/dsm_{placename}_cycle_{cycle}.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Make an animated GIF of changing ice surface from the PNG files\n",
    "# !convert -delay 120 -loop 0 figures/{placename}/dsm_*.png {gif_fname}\n",
    "gif_fname: str = (\n",
    "    f\"figures/{placename}/dsm_{placename}_cycles_{cycles[0]}-{cycles[-1]}.gif\"\n",
    ")\n",
    "subprocess.check_call(\n",
    "    [\n",
    "        \"convert\",\n",
    "        \"-delay\",\n",
    "        \"120\",\n",
    "        \"-loop\",\n",
    "        \"0\",\n",
    "        f\"figures/{placename}/dsm_*cycle_*.png\",\n",
    "        gif_fname,\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# HvPlot 2D interactive view of ice surface elevation grids over each ICESat-2 cycle\n",
    "dashboard: pn.layout.Column = pn.Column(\n",
    "    ds_lake.hvplot.image(x=\"x\", y=\"y\", clim=z_limits, cmap=\"gist_earth\", data_aspect=1)\n",
    "    # * ds_lake.hvplot.contour(x=\"x\", y=\"y\", clim=z_limits, data_aspect=1)\n",
    ")\n",
    "dashboard.show(port=30227)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Along track plots of ice surface elevation change over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Select a few Reference Ground tracks to look at\n",
    "savefig_tasks: list = []  # empty list of save figure tasks\n",
    "rgts: list = [int(rgt) for rgt in lake.refgtracks.split(\"|\")]\n",
    "print(f\"Looking at Reference Ground Tracks: {rgts}\")\n",
    "track_dict: dict = deepicedrain.split_tracks(df=df_lake.to_pandas())\n",
    "for rgtpair, df_ in track_dict.items():\n",
    "    # Transect plot along a reference ground track\n",
    "    fig = dask.delayed(obj=deepicedrain.plot_alongtrack)(\n",
    "        df=df_, rgtpair=rgtpair, regionname=region.name, oldtonew=draining\n",
    "    )\n",
    "    savefig_task = fig.savefig(\n",
    "        fname=f\"figures/{placename}/alongtrack_{placename}_{rgtpair}.png\"\n",
    "    )\n",
    "    savefig_tasks.append(savefig_task)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "futures = [client.compute(savefig_task) for savefig_task in savefig_tasks]\n",
    "for _ in tqdm.tqdm(\n",
    "    iterable=dask.distributed.as_completed(futures=futures), total=len(savefig_tasks)\n",
    "):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "encoding": "# -*- coding: utf-8 -*-",
   "formats": "ipynb,py:hydrogen"
  },
  "kernelspec": {
   "display_name": "deepicedrain",
   "language": "python",
   "name": "deepicedrain"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}