{
"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
}