atlxi_xover.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "# **ICESat-2 Crossover Track Analysis**\n",
    "\n",
    "To increase the temporal resolution of\n",
    "our ice elevation change analysis\n",
    "(i.e. at time periods less than\n",
    "the 91 day repeat cycle of ICESat-2),\n",
    "we can look at the locations where the\n",
    "ICESat-2 tracks intersect and get the\n",
    "height values there!\n",
    "Uses [pygmt.x2sys_cross](https://www.pygmt.org/v0.2.0/api/generated/pygmt.x2sys_cross.html).\n",
    "\n",
    "References:\n",
    "- Wessel, P. (2010). Tools for analyzing intersecting tracks: The x2sys package.\n",
    "Computers & Geosciences, 36(3), 348–354. https://doi.org/10.1016/j.cageo.2009.05.009"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import os\n",
    "\n",
    "import dask\n",
    "import deepicedrain\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pint\n",
    "import pint_pandas\n",
    "import pygmt\n",
    "import shapely.geometry\n",
    "import tqdm\n",
    "import uncertainties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "ureg = pint.UnitRegistry()\n",
    "pint_pandas.PintType.ureg = ureg\n",
    "tag: str = \"X2SYS\"\n",
    "os.environ[\"X2SYS_HOME\"] = os.path.abspath(tag)\n",
    "client = dask.distributed.Client(\n",
    "    n_workers=4, threads_per_worker=1, env={\"X2SYS_HOME\": os.environ[\"X2SYS_HOME\"]}\n",
    ")\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_date, max_date = (\"2019-03-29\", \"2020-12-24\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize X2SYS database in the X2SYS/ICESAT2 folder\n",
    "pygmt.x2sys_init(\n",
    "    tag=\"ICESAT2\",\n",
    "    fmtfile=f\"{tag}/ICESAT2/xyht\",\n",
    "    suffix=\"tsv\",\n",
    "    units=[\"de\", \"se\"],  # distance in metres, speed in metres per second\n",
    "    gap=\"d250e\",  # distance gap up to 250 metres allowed\n",
    "    force=True,\n",
    "    verbose=\"q\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Select a subglacial lake to examine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Save or load dhdt data from Parquet file\n",
    "placename: str = \"whillans_upstream\"  # \"slessor_downstream\"\n",
    "df_dhdt: pd.DataFrame = pd.read_parquet(f\"ATLXI/df_dhdt_{placename.lower()}.parquet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Subset data to lake of interest\n",
    "placename: str = region.name.lower().replace(\" \", \"_\")\n",
    "df_lake: pd.DataFrame = region.subset(data=df_dhdt)  # bbox subset\n",
    "gdf_lake = gpd.GeoDataFrame(\n",
    "    df_lake, geometry=gpd.points_from_xy(x=df_lake.x, y=df_lake.y, crs=3031)\n",
    ")\n",
    "df_lake: pd.DataFrame = df_lake.loc[gdf_lake.within(lake.geometry)]  # polygon subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Run crossover analysis on all tracks\n",
    "track_dict: dict = deepicedrain.split_tracks(df=df_lake)\n",
    "rgts, tracks = track_dict.keys(), track_dict.values()\n",
    "# Parallelized paired crossover analysis\n",
    "futures: list = []\n",
    "for rgt1, rgt2 in itertools.combinations(rgts, r=2):\n",
    "    # skip if same referencegroundtrack but different laser pair\n",
    "    # as they are parallel and won't cross\n",
    "    if rgt1[:4] == rgt2[:4]:\n",
    "        continue\n",
    "    track1 = track_dict[rgt1][[\"x\", \"y\", \"h_corr\", \"utc_time\"]]\n",
    "    track2 = track_dict[rgt2][[\"x\", \"y\", \"h_corr\", \"utc_time\"]]\n",
    "    shape1 = shapely.geometry.LineString(coordinates=track1[[\"x\", \"y\"]].to_numpy())\n",
    "    shape2 = shapely.geometry.LineString(coordinates=track2[[\"x\", \"y\"]].to_numpy())\n",
    "    if not shape1.intersects(shape2):\n",
    "        continue\n",
    "    future = client.submit(\n",
    "        key=f\"{rgt1}x{rgt2}\",\n",
    "        func=pygmt.x2sys_cross,\n",
    "        tracks=[track1, track2],\n",
    "        tag=\"ICESAT2\",\n",
    "        # region=[-460000, -400000, -560000, -500000],\n",
    "        interpolation=\"l\",  # linear interpolation\n",
    "        coe=\"e\",  # external crossovers\n",
    "        trackvalues=True,  # Get track 1 height (h_1) and track 2 height (h_2)\n",
    "        # trackvalues=False,  # Get crossover error (h_X) and mean height value (h_M)\n",
    "        # outfile=\"xover_236_562.tsv\"\n",
    "    )\n",
    "    futures.append(future)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "crossovers: dict = {}\n",
    "for f in tqdm.tqdm(\n",
    "    iterable=dask.distributed.as_completed(futures=futures), total=len(futures)\n",
    "):\n",
    "    if f.status != \"error\":  # skip those track pairs which don't intersect\n",
    "        crossovers[f.key] = f.result().dropna().reset_index(drop=True)\n",
    "\n",
    "df_cross: pd.DataFrame = pd.concat(objs=crossovers, names=[\"track1_track2\", \"id\"])\n",
    "df: pd.DataFrame = df_cross.reset_index(level=\"track1_track2\").reset_index(drop=True)\n",
    "# Report on how many unique crossover intersections there were\n",
    "# df.plot.scatter(x=\"x\", y=\"y\")  # quick plot of our crossover points\n",
    "print(\n",
    "    f\"{len(df.groupby(by=['x', 'y']))} crossover intersection point locations found \"\n",
    "    f\"with {len(df)} crossover height-time pairs \"\n",
    "    f\"over {len(tracks)} tracks\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate crossover error\n",
    "df[\"h_X\"]: pd.Series = df.h_2 - df.h_1  # crossover error (i.e. height difference)\n",
    "df[\"t_D\"]: pd.Series = df.t_2 - df.t_1  # elapsed time in ns (i.e. time difference)\n",
    "ns_in_yr: int = 365.25 * 24 * 60 * 60 * 1_000_000_000  # nanoseconds in a year\n",
    "df[\"dhdt\"]: pd.Series = df.h_X / (df.t_D.astype(np.int64) / ns_in_yr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Get some summary statistics of our crossover errors\n",
    "sumstats: pd.DataFrame = df[[\"h_X\", \"t_D\", \"dhdt\"]].describe()\n",
    "# Find location with highest absolute crossover error, and most sudden height change\n",
    "max_h_X: pd.Series = df.iloc[np.nanargmax(df.h_X.abs())]  # highest crossover error\n",
    "max_dhdt: pd.Series = df.iloc[df.dhdt.argmax()]  # most sudden change in height"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2D Map view of crossover points\n",
    "\n",
    "Bird's eye view of the crossover points\n",
    "overlaid on top of the ICESat-2 tracks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2D plot of crossover locations\n",
    "var: str = \"h_X\"\n",
    "fig = pygmt.Figure()\n",
    "# Setup basemap\n",
    "plotregion = pygmt.info(table=df[[\"x\", \"y\"]], spacing=1000)\n",
    "pygmt.makecpt(cmap=\"batlow\", series=[sumstats[var][\"25%\"], sumstats[var][\"75%\"]])\n",
    "# Map frame in metre units\n",
    "fig.basemap(frame=\"f\", region=plotregion, projection=\"X8c\")\n",
    "# Plot actual track points in green\n",
    "for track in tracks:\n",
    "    tracklabel = f\"{track.iloc[0].referencegroundtrack} {track.iloc[0].pairtrack}\"\n",
    "    fig.plot(\n",
    "        x=track.x,\n",
    "        y=track.y,\n",
    "        pen=\"thinnest,green,.\",\n",
    "        style=f'qN+1:+l\"{tracklabel}\"+f3p,Helvetica,darkgreen',\n",
    "    )\n",
    "# Plot crossover point locations\n",
    "fig.plot(x=df.x, y=df.y, color=df.h_X, cmap=True, style=\"c0.1c\", pen=\"thinnest\")\n",
    "# Plot lake boundary in blue\n",
    "lakex, lakey = lake.geometry.exterior.coords.xy\n",
    "fig.plot(x=lakex, y=lakey, pen=\"thin,blue,-.\")\n",
    "# Map frame in kilometre units\n",
    "fig.basemap(\n",
    "    frame=[\n",
    "        f'WSne+t\"Crossover points at {region.name}\"',\n",
    "        'xaf+l\"Polar Stereographic X (km)\"',\n",
    "        'yaf+l\"Polar Stereographic Y (km)\"',\n",
    "    ],\n",
    "    region=plotregion / 1000,\n",
    "    projection=\"X8c\",\n",
    ")\n",
    "fig.colorbar(position=\"JMR+e\", frame=['x+l\"Crossover Error\"', \"y+lm\"])\n",
    "fig.savefig(f\"figures/{placename}/crossover_area_{placename}_{min_date}_{max_date}.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot Crossover Elevation time-series\n",
    "\n",
    "Plot elevation change over time at:\n",
    "\n",
    "1. One single crossover point location\n",
    "2. Many crossover locations over an area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tidy up dataframe first using pd.wide_to_long\n",
    "# I.e. convert 't_1', 't_2', 'h_1', 'h_2' columns into just 't' and 'h'.\n",
    "df_th: pd.DataFrame = deepicedrain.wide_to_long(\n",
    "    df=df.loc[:, [\"track1_track2\", \"x\", \"y\", \"t_1\", \"t_2\", \"h_1\", \"h_2\"]],\n",
    "    stubnames=[\"t\", \"h\"],\n",
    "    j=\"track\",\n",
    ")\n",
    "df_th: pd.DataFrame = df_th.drop_duplicates(ignore_index=True)\n",
    "df_th: pd.DataFrame = df_th.sort_values(by=\"t\").reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot at single location with **maximum** absolute crossover height error (max_h_X)\n",
    "df_max = df_th.query(expr=\"x == @max_h_X.x & y == @max_h_X.y\")\n",
    "track1, track2 = df_max.track1_track2.iloc[0].split(\"x\")\n",
    "print(f\"{max_h_X.h_X:.2f} metres height change at {max_h_X.x}, {max_h_X.y}\")\n",
    "plotregion = np.array(\n",
    "    [df_max.t.min(), df_max.t.max(), *pygmt.info(table=df_max[[\"h\"]], spacing=2.5)[:2]]\n",
    ")\n",
    "plotregion += np.array([-pd.Timedelta(2, unit=\"W\"), +pd.Timedelta(2, unit=\"W\"), 0, 0])\n",
    "\n",
    "fig = pygmt.Figure()\n",
    "with pygmt.config(\n",
    "    FONT_ANNOT_PRIMARY=\"9p\", FORMAT_TIME_PRIMARY_MAP=\"abbreviated\", FORMAT_DATE_MAP=\"o\"\n",
    "):\n",
    "    fig.basemap(\n",
    "        projection=\"X12c/8c\",\n",
    "        region=plotregion,\n",
    "        frame=[\n",
    "            f'WSne+t\"Max elevation change over time at {region.name}\"',\n",
    "            \"pxa1Of1o+lDate\",  # primary time axis, 1 mOnth annotation and minor axis\n",
    "            \"sx1Y\",  # secondary time axis, 1 Year intervals\n",
    "            'yaf+l\"Elevation at crossover (m)\"',\n",
    "        ],\n",
    "    )\n",
    "fig.text(\n",
    "    text=f\"Track {track1} and {track2} crossover\",\n",
    "    position=\"TC\",\n",
    "    offset=\"jTC0c/0.2c\",\n",
    "    verbose=\"q\",\n",
    ")\n",
    "# Plot data points\n",
    "fig.plot(x=df_max.t, y=df_max.h, style=\"c0.15c\", color=\"darkblue\", pen=\"thin\")\n",
    "# Plot dashed line connecting points\n",
    "fig.plot(x=df_max.t, y=df_max.h, pen=f\"faint,blue,-\")\n",
    "fig.savefig(\n",
    "    f\"figures/{placename}/crossover_point_{placename}_{track1}_{track2}_{min_date}_{max_date}.png\"\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all crossover height points over time over the lake area\n",
    "fig = deepicedrain.plot_crossovers(df=df_th, regionname=region.name)\n",
    "fig.savefig(f\"figures/{placename}/crossover_many_{placename}_{min_date}_{max_date}.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate height anomaly at crossover point as\n",
    "# height at t=n minus height at t=0 (first observation date at crossover point)\n",
    "anomfunc = lambda h: h - h.iloc[0]  # lambda h: h - h.mean()\n",
    "df_th[\"h_anom\"] = df_th.groupby(by=\"track1_track2\").h.transform(func=anomfunc)\n",
    "# Calculate ice volume displacement (dvol) in unit metres^3\n",
    "# and rolling mean height anomaly (h_roll) in unit metres\n",
    "surface_area: pint.Quantity = lake.geometry.area * ureg.metre ** 2\n",
    "ice_dvol: pd.Series = deepicedrain.ice_volume_over_time(\n",
    "    df_elev=df_th.astype(dtype={\"h_anom\": \"pint[metre]\"}),\n",
    "    surface_area=surface_area,\n",
    "    time_col=\"t\",\n",
    "    outfile=f\"figures/{placename}/ice_dvol_dt_{placename}.txt\",\n",
    ")\n",
    "df_th[\"h_roll\"]: pd.Series = uncertainties.unumpy.nominal_values(\n",
    "    arr=ice_dvol.pint.magnitude / surface_area.magnitude\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all crossover height point anomalies over time over the lake area\n",
    "fig = deepicedrain.plot_crossovers(\n",
    "    df=df_th,\n",
    "    regionname=region.name,\n",
    "    elev_var=\"h_anom\",\n",
    "    outline_points=f\"figures/{placename}/{placename}.gmt\",\n",
    ")\n",
    "fig.plot(x=df_th.t, y=df_th.h_roll, pen=\"thick,-\")  # plot rolling mean height anomaly\n",
    "fig.savefig(\n",
    "    f\"figures/{placename}/crossover_anomaly_{placename}_{min_date}_{max_date}.png\"\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combined ice volume displacement plot\n",
    "\n",
    "Showing how subglacial water cascades down a drainage basin!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "fig.basemap(\n",
    "    region=f\"2019-02-28/2020-09-30/-0.3/0.5\",\n",
    "    frame=[\"wSnE\", \"xaf\", 'yaf+l\"Ice Volume Displacement (km@+3@+)\"'],\n",
    ")\n",
    "pygmt.makecpt(cmap=\"davosS\", color_model=\"+c\", series=(-2, 4, 0.5))\n",
    "for i, (_placename, linestyle) in enumerate(\n",
    "    iterable=zip(\n",
    "        [\"whillans_ix\", \"subglacial_lake_whillans\", \"lake_12\", \"whillans_7\"],\n",
    "        [\"\", \".-\", \"-\", \"..-\"],\n",
    "    )\n",
    "):\n",
    "    fig.plot(\n",
    "        data=f\"figures/{_placename}/ice_dvol_dt_{_placename}.txt\",\n",
    "        cmap=True,\n",
    "        pen=f\"thick,{linestyle}\",\n",
    "        zvalue=i,\n",
    "        label=_placename,\n",
    "        columns=\"0,3\",  # time column (0), ice_dvol column (3)\n",
    "    )\n",
    "fig.text(\n",
    "    position=\"TL\",\n",
    "    offset=\"j0.2c\",\n",
    "    text=\"Whillans Ice Stream Central Catchment active subglacial lakes\",\n",
    ")\n",
    "fig.legend(position=\"jML+jML+o0.2c\", box=\"+gwhite\")\n",
    "fig.savefig(\"figures/cascade_whillans_ice_stream.png\")\n",
    "fig.show()"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}