Skip to content

Commit

Permalink
💥 ICESat-2 elevation change at crossovers within a subglacial lake
Browse files Browse the repository at this point in the history 
Finally, showing (nicely) how sudden active subglacial lakes can fill up (~14m over 4 months at Slessor 23)! Made a plot_crossovers function that can plot both raw and normalized height changes over time, for every crossover location from intersecting ICESat-2 tracks. Refactored a lot of the code to make use of PyGMT v0.2.0 capabilities (e.g. nicer `pygmt.info` and `pygmt.x2sys_cross`), and also some pandas groupby magic. Also improved the quality of the vizplots.py docstring by including a fancy ASCII text representation of the expected input table data and output figure!
  • Loading branch information
@weiji14
weiji14 committed Sep 16, 2020
1 parent fc1eba2 commit d02726c
Show file tree
Hide file tree
Showing 6 changed files with 324 additions and 155 deletions.
147 changes: 79 additions & 68 deletions atlxi_lake.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
},
"outputs": [],
"source": [
"import itertools\n",
"import os\n",
"\n",
"os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"1\"\n",
Expand All @@ -53,13 +54,36 @@
"import deepicedrain"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tag: str = \"X2SYS\"\n",
"os.environ[\"X2SYS_HOME\"] = os.path.abspath(tag)\n",
"client = dask.distributed.Client(\n",
" n_workers=64, threads_per_worker=1, env={\"X2SYS_HOME\": os.environ[\"X2SYS_HOME\"]}\n",
")\n",
"client"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Data Preparation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"min_date, max_date = (\"2018-10-14\", \"2020-05-13\")"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand Down Expand Up @@ -271,7 +295,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"# Plot clusters on a map in colour, noise points/outliers as small dots\n",
Expand Down Expand Up @@ -349,7 +375,7 @@
"outputs": [],
"source": [
"# Choose one draining/filling lake\n",
"draining: bool = True # False\n",
"draining: bool = False # False\n",
"placename: str = \"Slessor\" # \"Whillans\"\n",
"lakes: gpd.GeoDataFrame = antarctic_lakes.query(expr=\"basin_name == @placename\")\n",
"lake = lakes.loc[lakes.maxabsdhdt.idxmin() if draining else lakes.maxabsdhdt.idxmax()]\n",
Expand Down Expand Up @@ -427,7 +453,7 @@
"we can look at the locations where the\n",
"ICESat-2 tracks intersect and get the\n",
"height values there!\n",
"Uses [x2sys_cross](https://docs.generic-mapping-tools.org/6.1/supplements/x2sys/x2sys_cross).\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",
Expand All @@ -441,9 +467,6 @@
"outputs": [],
"source": [
"# Initialize X2SYS database in the X2SYS/ICESAT2 folder\n",
"tag = \"X2SYS\"\n",
"os.environ[\"X2SYS_HOME\"] = os.path.abspath(tag)\n",
"os.getcwd()\n",
"pygmt.x2sys_init(\n",
" tag=\"ICESAT2\",\n",
" fmtfile=f\"{tag}/ICESAT2/xyht\",\n",
Expand All @@ -464,20 +487,18 @@
"outputs": [],
"source": [
"# Run crossover analysis on all tracks\n",
"rgts: list = [135, 327, 388, 577, 1080, 1272] # Whillans upstream\n",
"# rgts: list = [236, 501, 562, 1181] # Whillans_downstream\n",
"tracks = [f\"{tag}/track_{i}.tsv\" for i in rgts]\n",
"assert all(os.path.exists(k) for k in tracks)\n",
"\n",
"rgts, tracks = track_dict.keys(), track_dict.values()\n",
"# Parallelized paired crossover analysis\n",
"futures: list = []\n",
"for track1, track2 in itertools.combinations(rgts, r=2):\n",
"for rgt1, rgt2 in itertools.combinations(rgts, r=2):\n",
" track1 = track_dict[rgt1][[\"x\", \"y\", \"h_corr\", \"utc_time\"]]\n",
" track2 = track_dict[rgt2][[\"x\", \"y\", \"h_corr\", \"utc_time\"]]\n",
" future = client.submit(\n",
" key=f\"{track1}_{track2}\",\n",
" key=f\"{rgt1}x{rgt2}\",\n",
" func=pygmt.x2sys_cross,\n",
" tracks=[f\"{tag}/track_{track1}.tsv\", f\"{tag}/track_{track2}.tsv\"],\n",
" tracks=[track1, track2],\n",
" tag=\"ICESAT2\",\n",
" region=[-460000, -400000, -560000, -500000],\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",
Expand Down Expand Up @@ -563,38 +584,40 @@
"var: str = \"h_X\"\n",
"fig = pygmt.Figure()\n",
"# Setup basemap\n",
"region = np.array([df.x.min(), df.x.max(), df.y.min(), df.y.max()])\n",
"buffer = np.array([-2000, +2000, -2000, +2000])\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=region + buffer, projection=\"X8c\")\n",
"fig.basemap(frame=\"f\", region=plotregion, projection=\"X8c\")\n",
"# Plot actual track points\n",
"for track in tracks:\n",
" fig.plot(data=track, color=\"green\", style=\"c0.01c\")\n",
" fig.plot(x=track.x, y=track.y, color=\"green\", style=\"c0.01c\")\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\n",
"lakex, lakey = lake.geometry.exterior.coords.xy\n",
"fig.plot(x=lakex, y=lakey, pen=\"thin\")\n",
"# Map frame in kilometre units\n",
"fig.basemap(\n",
" frame=[\n",
" \"WSne\",\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=(region + buffer) / 1000,\n",
" region=plotregion / 1000,\n",
" projection=\"X8c\",\n",
")\n",
"fig.colorbar(position=\"JMR\", frame=['x+l\"Crossover Error\"', \"y+lm\"])\n",
"fig.savefig(\"figures/crossover_area.png\")\n",
"fig.savefig(f\"figures/crossover_area_{placename}.png\")\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1D plots of height changing over time\n",
"### Plot Crossover Elevation time-series\n",
"\n",
"Plot height change over time at:\n",
"Plot elevation change over time at:\n",
"\n",
"1. One single crossover point location\n",
"2. Many crossover locations over an area"
Expand All @@ -609,7 +632,7 @@
"# 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[[\"track1_track2\", \"x\", \"y\", \"t_1\", \"t_2\", \"h_1\", \"h_2\"]],\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",
Expand All @@ -622,24 +645,24 @@
"metadata": {},
"outputs": [],
"source": [
"# 1D Plot at location with **maximum** absolute crossover height error (max_h_X)\n",
"# 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\").sort_values(by=\"t\")\n",
"track1, track2 = df_max.track1_track2.iloc[0].split(\"_\")\n",
"track1, track2 = df_max.track1_track2.iloc[0].split(\"x\")\n",
"print(f\"{round(max_h_X.h_X, 2)} metres height change at {max_h_X.x}, {max_h_X.y}\")\n",
"t_min = (df_max.t.min() - pd.Timedelta(2, unit=\"W\")).isoformat()\n",
"t_max = (df_max.t.max() + pd.Timedelta(2, unit=\"W\")).isoformat()\n",
"h_min = df_max.h.min() - 0.2\n",
"h_max = df_max.h.max() + 0.4\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=[t_min, t_max, h_min, h_max],\n",
" region=plotregion,\n",
" frame=[\n",
" \"WSne\",\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",
Expand All @@ -655,7 +678,9 @@
"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(f\"figures/crossover_{track1}_{track2}_{min_date}_{max_date}.png\")\n",
"fig.savefig(\n",
" f\"figures/crossover_point_{placename}_{track1}_{track2}_{min_date}_{max_date}.png\"\n",
")\n",
"fig.show()"
]
},
Expand All @@ -665,41 +690,27 @@
"metadata": {},
"outputs": [],
"source": [
"# 1D plots of a crossover area, all the height points over time\n",
"t_min = (df_th.t.min() - pd.Timedelta(1, unit=\"W\")).isoformat()\n",
"t_max = (df_th.t.max() + pd.Timedelta(1, unit=\"W\")).isoformat()\n",
"h_min = df_th.h.min() - 0.2\n",
"h_max = df_th.h.max() + 0.2\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/12c\",\n",
" region=[t_min, t_max, h_min, h_max],\n",
" frame=[\n",
" \"WSne\",\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",
"# Plot all crossover height points over time over the lake area\n",
"fig = deepicedrain.vizplots.plot_crossovers(df=df_th, regionname=region.name)\n",
"fig.savefig(f\"figures/crossover_many_{placename}_{min_date}_{max_date}.png\")\n",
"fig.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot all crossover height points over time over the lake area\n",
"# with height values normalized to 0 from the first observation date\n",
"normfunc = lambda h: h - h.iloc[0] # lambda h: h - h.mean()\n",
"df_th[\"h_norm\"] = df_th.groupby(by=\"track1_track2\").h.transform(func=normfunc)\n",
"\n",
"crossovers = df_th.groupby(by=[\"x\", \"y\"])\n",
"pygmt.makecpt(cmap=\"categorical\", series=[1, len(crossovers) + 1, 1])\n",
"for i, ((x_coord, y_coord), indexes) in enumerate(crossovers.indices.items()):\n",
" df_ = df_th.loc[indexes].sort_values(by=\"t\")\n",
" # if df_.h.max() - df_.h.min() > 1.0: # plot only > 1 metre height change\n",
" track1, track2 = df_.track1_track2.iloc[0].split(\"_\")\n",
" label = f'\"Track {track1} {track2}\"'\n",
" fig.plot(x=df_.t, y=df_.h, Z=i, style=\"c0.1c\", cmap=True, pen=\"thin+z\", label=label)\n",
" # Plot line connecting points\n",
" fig.plot(\n",
" x=df_.t, y=df_.h, Z=i, pen=f\"faint,+z,-\", cmap=True\n",
" ) # , label=f'\"+g-1l+s0.15c\"')\n",
"fig.legend(position=\"JMR+JMR+o0.2c\", box=\"+gwhite+p1p\")\n",
"fig.savefig(f\"figures/crossover_many_{min_date}_{max_date}.png\")\n",
"fig = deepicedrain.vizplots.plot_crossovers(\n",
" df=df_th, regionname=region.name, elev_var=\"h_norm\"\n",
")\n",
"fig.savefig(f\"figures/crossover_many_normalized_{placename}_{min_date}_{max_date}.png\")\n",
"fig.show()"
]
},
Expand Down
Loading

0 comments on commit d02726c

Please sign in to comment.