From 69e5b634544dbb58531bba75102034eca91b1861 Mon Sep 17 00:00:00 2001 From: Wei Ji Date: Wed, 14 Oct 2020 15:56:18 +1300 Subject: [PATCH 1/5] :sparkles: Spatiotemporal data cube of elev surfaces per ICESat-2 cycle It's easier to look at raster grid surfaces than point clouds, and visualize them in 3D too! Using the good ol' `blockmedian` and `surface` modules from PyGMT to interpolate ICESat-2 points onto a grid. This is done on a (3 month) cycle by cycle basis, and the resulting NetCDF files are stacked into a single NetCDF data cube using xarray. Will need to work on using a proper datetime64 time axis instead of just cycle number, and perhaps allow for kwargs into the `surface` instead of the current hardcoded parameters. --- deepicedrain/README.md | 1 + deepicedrain/__init__.py | 1 + deepicedrain/spatiotemporal.py | 119 +++++++++++++++++- .../tests/test_spatiotemporal_cube.py | 61 +++++++++ 4 files changed, 180 insertions(+), 2 deletions(-) create mode 100644 deepicedrain/tests/test_spatiotemporal_cube.py diff --git a/deepicedrain/README.md b/deepicedrain/README.md index 63a2e7f..17c6af7 100644 --- a/deepicedrain/README.md +++ b/deepicedrain/README.md @@ -17,6 +17,7 @@ Contents: - Region - Bounding box data class structure that has xarray subsetting capabilities and more! - deltatime_to_utctime - Converts GPS time from an epoch (default is 2018 Jan 1st) to UTC time - lonlat_to_xy - Reprojects longitude/latitude EPSG:4326 coordinates to x/y EPSG:3031 coordinates + - spatiotemporal_cube - Interpolates a time-series point cloud into an xarray.Dataset data cube - :card_file_box: extraload.py - Convenience functions for extracting, transforming and loading data - array_to_dataframe - Turns a 1D/2D numpy/dask array into a tidy pandas/dask dataframe table diff --git a/deepicedrain/__init__.py b/deepicedrain/__init__.py index 7b40fc2..2a68768 100644 --- a/deepicedrain/__init__.py +++ b/deepicedrain/__init__.py @@ -11,6 +11,7 @@ deltatime_to_utctime, lonlat_to_xy, point_in_polygon_gpu, + spatiotemporal_cube, ) from deepicedrain.vizplots import IceSat2Explorer, plot_alongtrack, plot_crossovers diff --git a/deepicedrain/spatiotemporal.py b/deepicedrain/spatiotemporal.py index 8ce5b86..191d158 100644 --- a/deepicedrain/spatiotemporal.py +++ b/deepicedrain/spatiotemporal.py @@ -5,16 +5,16 @@ import dataclasses import datetime import os +import shutil import tempfile +import datashader import geopandas as gpd import numpy as np import pandas as pd import pyproj import xarray as xr -import datashader - @dataclasses.dataclass(frozen=True) class Region: @@ -298,3 +298,118 @@ def point_in_polygon_gpu( point_labels.loc[poly_labels[label]] = poly_df_.loc[label][poly_label_col] return point_labels + + +def spatiotemporal_cube( + table: pd.DataFrame, + placename: str = "", + x_var: str = "x", + y_var: str = "y", + z_var: str = "h_corr", + spacing: int = 250, + cycles: list = None, + folder: str = "", +) -> xr.Dataset: + """ + Interpolates a time-series point cloud into an xarray.Dataset data cube. + Uses `pygmt`'s blockmedian and surface algorithms to produce individual + NetCDF grids, and `xarray` to stack each NetCDF grid into one dataset. + + Steps are as follows: + + 1. Create several xarray.DataArray grid surfaces from a table of points, + one for each time cycle. + 2. Stacked the grids along a time cycle axis into a xarray.Dataset which is + a spatiotemporal data cube with 'x', 'y' and 'cycle_number' dimensions. + + _1__2__3_ + * * / / / /| + * * / / / / | + * * * /__/__/__/ | y + * * * --> | | | | | + * * * | | | | / + * * |__|__|__|/ x + cycle + + Parameters + ---------- + table : pandas.DataFrame + A table containing the ICESat-2 track data from multiple cycles. It + should ideally have geographical columns called 'x', 'y', and attribute + columns like 'h_corr_1', 'h_corr_2', etc for each cycle time. + placename : str + Optional. A descriptive placename for the data (e.g. some_ice_stream), + to be used in the temporary NetCDF filename. + x_var : str + The x coordinate column name to use from the table data. Default is + 'x'. + y_var : str + The y coordinate column name to use from the table data. Default is + 'y'. + z_var : str + The z column name to use from the table data. This will be the + attribute that the surface algorithm will run on. Default is 'h_corr'. + spacing : float or str + The spatial resolution of the resulting grid, provided as a number or + as 'dx/dy' increments. This is passed on to `pygmt.blockmedian` and + `pygmt.surface`. Default is 250 (metres). + cycles : list + The cycle numbers to run the gridding algorithm on, e.g. [3, 4] will + use columns 'h_corr_3' and 'h_corr_4'. Default is None which will + automatically determine the cycles for a given z_var. + folder : str + The folder to keep the intermediate NetCDF file in. Default is to place + the files in the current working directory. + + Returns + ------- + cube : xarray.Dataset + A 3-dimensional data cube made of digital surfaces stacked along a time + cycle axis. + + """ + import pygmt + import tqdm + + # Determine grid's bounding box region (xmin, xmax, ymin, ymax) + grid_region: np.ndarray = pygmt.info( + table=table[[x_var, y_var]], spacing=f"s{spacing}" + ) + + # Create one grid surface for each time cycle + if cycles is None: + cycles: list = [ + int(col[len(z_var) + 1 :]) for col in table.columns if col.startswith(z_var) + ] + _placename = f"_{placename}" if placename else "" + for cycle in tqdm.tqdm(iterable=cycles): + df_trimmed = pygmt.blockmedian( + table=table[[x_var, y_var, f"{z_var}_{cycle}"]].dropna(), + region=grid_region, + spacing=f"{spacing}+e", + ) + outfile = f"{z_var}{_placename}_cycle_{cycle}.nc" + pygmt.surface( + data=df_trimmed.values, + region=grid_region, + spacing=spacing, + T=0.35, # tension factor + V="e", # error messages only + outfile=outfile, + ) + # print(pygmt.grdinfo(outfile)) + + # Stack several NetCDF grids into one NetCDF along the time cycle axis + paths: list = [f"{z_var}{_placename}_cycle_{cycle}.nc" for cycle in cycles] + dataset: xr.Dataset = xr.open_mfdataset( + paths=paths, + combine="nested", + concat_dim=[pd.Index(data=cycles, name="cycle_number")], + attrs_file=paths[-1], + ) + + # Move files into new folder if requested + if folder: + [shutil.move(src=path, dst=folder) for path in paths] + + return dataset diff --git a/deepicedrain/tests/test_spatiotemporal_cube.py b/deepicedrain/tests/test_spatiotemporal_cube.py new file mode 100644 index 0000000..84ff13a --- /dev/null +++ b/deepicedrain/tests/test_spatiotemporal_cube.py @@ -0,0 +1,61 @@ +""" +Tests creation of spatiotemporal NetCDF data cubes +""" +import os +import tempfile + +import pandas as pd +import pytest +import xarray as xr + +from deepicedrain import catalog, ndarray_to_parquet, spatiotemporal_cube + + +@pytest.fixture(scope="module", name="table") +def fixture_table(): + """ + Loads the sample ICESat-2 ATL11 data, and processes it into an suitable + pandas.DataFrame table format. + """ + dataset: xr.Dataset = catalog.test_data.atl11_test_case.to_dask() + with tempfile.TemporaryDirectory() as tmpdir: + parquetpath: str = os.path.join(tmpdir, "temp.parquet") + table: pd.DataFrame = ndarray_to_parquet( + ndarray=dataset, + parquetpath=parquetpath, + variables=["longitude", "latitude", "h_corr"], + ) + return table + + +def test_spatiotemporal_cube(table): + """ + Tests that creating a spatiotemporal NetCDF data cube works + """ + grid: xr.Dataset = spatiotemporal_cube( + table=table, + placename="greenland", + x_var="longitude", + y_var="latitude", + spacing=0.1, + folder=tempfile.gettempdir(), + ) + + assert isinstance(grid, xr.Dataset) + assert grid.dims == {"x": 8, "y": 16, "cycle_number": 2} + xr.testing.assert_allclose(a=grid.z.min(), b=xr.DataArray(data=-19568.74)) + xr.testing.assert_allclose(a=grid.z.max(), b=xr.DataArray(data=21167.637)) + xr.testing.assert_allclose( + a=grid.z.median(axis=(1, 2)), + b=xr.DataArray( + data=[764.9603, 764.9551], coords=[[1, 2]], dims=["cycle_number"] + ), + ) + assert "-Gh_corr_greenland_cycle_2.nc" in grid.history + + paths = [ + os.path.join(tempfile.gettempdir(), f"h_corr_greenland_cycle_{cycle}.nc") + for cycle in (1, 2) + ] + assert all(os.path.exists(path=path) for path in paths) + [os.remove(path=path) for path in paths] From ec1eac171962eb68ff85269a15a383f89b9792d5 Mon Sep 17 00:00:00 2001 From: Wei Ji Date: Wed, 21 Oct 2020 15:51:27 +1300 Subject: [PATCH 2/5] :recycle: Better handling of file move and proj4 in spatiotemporal_cube Refactor to move the NetCDF files into the folder correctly, overwriting the filepath(s) if the file is already present. Also added a 'projection' parameter so that the projection information is stored in the NetCDF grid and viewable directly in QGIS without needing to set a projection first. Using the proj4 string of Antarctic Polar Stereographic EPSG:3031 from https://spatialreference.org/ref/epsg/3031/proj4/. Note to self, the grid_mapping coordinate should actually be an attribute, either handle that using xarray.open_rasterio, use rioxarray, or wait until the pangeo/xarray community standardizes on what to do with this. --- deepicedrain/spatiotemporal.py | 19 ++++++++++++++----- .../tests/test_spatiotemporal_cube.py | 2 +- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/deepicedrain/spatiotemporal.py b/deepicedrain/spatiotemporal.py index 191d158..3bbbbeb 100644 --- a/deepicedrain/spatiotemporal.py +++ b/deepicedrain/spatiotemporal.py @@ -308,6 +308,7 @@ def spatiotemporal_cube( z_var: str = "h_corr", spacing: int = 250, cycles: list = None, + projection: str = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", folder: str = "", ) -> xr.Dataset: """ @@ -357,6 +358,11 @@ def spatiotemporal_cube( The cycle numbers to run the gridding algorithm on, e.g. [3, 4] will use columns 'h_corr_3' and 'h_corr_4'. Default is None which will automatically determine the cycles for a given z_var. + projection : str + The proj4 string to store in the NetCDF output, will be passed directly + to `pygmt.surface`'s J (projection) argument. Default is '+proj=stere + +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 + +units=m +no_defs', i.e. Antarctic Polar Stereographic EPSG:3031. folder : str The folder to keep the intermediate NetCDF file in. Default is to place the files in the current working directory. @@ -393,14 +399,21 @@ def spatiotemporal_cube( data=df_trimmed.values, region=grid_region, spacing=spacing, + J=f'"{projection}"', # projection T=0.35, # tension factor V="e", # error messages only outfile=outfile, ) # print(pygmt.grdinfo(outfile)) - # Stack several NetCDF grids into one NetCDF along the time cycle axis + # Move files into new folder if requested paths: list = [f"{z_var}{_placename}_cycle_{cycle}.nc" for cycle in cycles] + if folder: + paths: list = [ + shutil.move(src=path, dst=os.path.join(folder, path)) for path in paths + ] + + # Stack several NetCDF grids into one NetCDF along the time cycle axis dataset: xr.Dataset = xr.open_mfdataset( paths=paths, combine="nested", @@ -408,8 +421,4 @@ def spatiotemporal_cube( attrs_file=paths[-1], ) - # Move files into new folder if requested - if folder: - [shutil.move(src=path, dst=folder) for path in paths] - return dataset diff --git a/deepicedrain/tests/test_spatiotemporal_cube.py b/deepicedrain/tests/test_spatiotemporal_cube.py index 84ff13a..7d8bd8d 100644 --- a/deepicedrain/tests/test_spatiotemporal_cube.py +++ b/deepicedrain/tests/test_spatiotemporal_cube.py @@ -42,7 +42,7 @@ def test_spatiotemporal_cube(table): ) assert isinstance(grid, xr.Dataset) - assert grid.dims == {"x": 8, "y": 16, "cycle_number": 2} + assert grid.dims == {"x": 8, "y": 16, "cycle_number": 2, "grid_mapping": 12} xr.testing.assert_allclose(a=grid.z.min(), b=xr.DataArray(data=-19568.74)) xr.testing.assert_allclose(a=grid.z.max(), b=xr.DataArray(data=21167.637)) xr.testing.assert_allclose( From 665995ad26a857bbb5e60938bd62fb7ceff61069 Mon Sep 17 00:00:00 2001 From: Wei Ji Date: Mon, 26 Oct 2020 13:23:06 +1300 Subject: [PATCH 3/5] :hammer: Impose 3 median absolute deviation limit on grid surface output MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Better handle elevation outliers that had resulted in sharp peaks and troughs in the output NetCDF grid. Doing so by telling `pygmt.surface` to limit the output solution to a lower and upper bound, determined as ±3 median absolute deviations from the median of the entire dataframe table. Originally I was clipping the grid in a post-processing step, but decided it's much better to impose the limits in the surface algorithm itself! The elevation statistics of the test case in Greenland also looks more sensible now than before. --- deepicedrain/spatiotemporal.py | 22 ++++++++++++++++++- .../tests/test_spatiotemporal_cube.py | 6 ++--- 2 files changed, 24 insertions(+), 4 deletions(-) diff --git a/deepicedrain/spatiotemporal.py b/deepicedrain/spatiotemporal.py index 3bbbbeb..ff52d65 100644 --- a/deepicedrain/spatiotemporal.py +++ b/deepicedrain/spatiotemporal.py @@ -13,6 +13,7 @@ import numpy as np import pandas as pd import pyproj +import scipy.stats import xarray as xr @@ -307,6 +308,7 @@ def spatiotemporal_cube( y_var: str = "y", z_var: str = "h_corr", spacing: int = 250, + clip_limits: bool = True, cycles: list = None, projection: str = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", folder: str = "", @@ -354,6 +356,10 @@ def spatiotemporal_cube( The spatial resolution of the resulting grid, provided as a number or as 'dx/dy' increments. This is passed on to `pygmt.blockmedian` and `pygmt.surface`. Default is 250 (metres). + clip_limits : bool + Whether or not to clip the output grid surface to ± 3 times the median + absolute deviation of the data table's z-values. Useful for handling + outlier values in the data table. Default is True (will clip). cycles : list The cycle numbers to run the gridding algorithm on, e.g. [3, 4] will use columns 'h_corr_3' and 'h_corr_4'. Default is None which will @@ -382,11 +388,24 @@ def spatiotemporal_cube( table=table[[x_var, y_var]], spacing=f"s{spacing}" ) - # Create one grid surface for each time cycle + # Automatically determine list of cycles if None is given if cycles is None: cycles: list = [ int(col[len(z_var) + 1 :]) for col in table.columns if col.startswith(z_var) ] + + # Limit surface output to within 3 median absolute deviations of median value + if clip_limits: + z_values = table[[f"{z_var}_{cycle}" for cycle in cycles]] + median: float = np.nanmedian(z_values) + meddev: float = scipy.stats.median_abs_deviation( + x=z_values, axis=None, nan_policy="omit" + ) + limits: list = [f"l{median - 3 * meddev}", f"u{median + 3 * meddev}"] + else: + limits = None + + # Create one grid surface for each time cycle _placename = f"_{placename}" if placename else "" for cycle in tqdm.tqdm(iterable=cycles): df_trimmed = pygmt.blockmedian( @@ -400,6 +419,7 @@ def spatiotemporal_cube( region=grid_region, spacing=spacing, J=f'"{projection}"', # projection + L=limits, # lower and upper limits T=0.35, # tension factor V="e", # error messages only outfile=outfile, diff --git a/deepicedrain/tests/test_spatiotemporal_cube.py b/deepicedrain/tests/test_spatiotemporal_cube.py index 7d8bd8d..d879fc8 100644 --- a/deepicedrain/tests/test_spatiotemporal_cube.py +++ b/deepicedrain/tests/test_spatiotemporal_cube.py @@ -43,12 +43,12 @@ def test_spatiotemporal_cube(table): assert isinstance(grid, xr.Dataset) assert grid.dims == {"x": 8, "y": 16, "cycle_number": 2, "grid_mapping": 12} - xr.testing.assert_allclose(a=grid.z.min(), b=xr.DataArray(data=-19568.74)) - xr.testing.assert_allclose(a=grid.z.max(), b=xr.DataArray(data=21167.637)) + xr.testing.assert_allclose(a=grid.z.min(), b=xr.DataArray(data=1435.1884)) + xr.testing.assert_allclose(a=grid.z.max(), b=xr.DataArray(data=1972.5968)) xr.testing.assert_allclose( a=grid.z.median(axis=(1, 2)), b=xr.DataArray( - data=[764.9603, 764.9551], coords=[[1, 2]], dims=["cycle_number"] + data=[1655.0094, 1654.4307], coords=[[1, 2]], dims=["cycle_number"] ), ) assert "-Gh_corr_greenland_cycle_2.nc" in grid.history From a63fad20c2eda57e8aaad95d53c0fbe353772204 Mon Sep 17 00:00:00 2001 From: Wei Ji Date: Mon, 26 Oct 2020 17:50:38 +1300 Subject: [PATCH 4/5] :triangular_flag_on_post: Refresh crossover track plotting scripts Main changes are an elev_filter setting in `deepicedrain.vizplots.plot_crossovers`, and the addition of annotated track numbers in the crossover_area map plot using GMT's quoted line functionality. Lots of other little bits too, including a refresh of the subglacial_lake ID (key) and name (value) pairs, a quick improvement to the inefficient brute-force crossover for-loop, and organizing figures inside 'figures/{lakename}' instead of just the 'figures' directory. Still a big chunk of local code yet to be committed, need to refactor and streamline the processing to be able to loop through many subglacial lakes in a more hands off manner! --- atlxi_lake.ipynb | 81 +++++++++++++++++++++++++++------------- atlxi_lake.py | 77 ++++++++++++++++++++++++++------------ deepicedrain/vizplots.py | 11 ++++-- 3 files changed, 116 insertions(+), 53 deletions(-) diff --git a/atlxi_lake.ipynb b/atlxi_lake.ipynb index fe1aa3a..2e8bdb6 100644 --- a/atlxi_lake.ipynb +++ b/atlxi_lake.ipynb @@ -917,11 +917,14 @@ "outputs": [], "source": [ "# Save or load dhdt data from Parquet file\n", - "placename: str = \"Recovery\" # \"Whillans\"\n", - "drainage_basins: gpd.GeoDataFrame = drainage_basins.set_index(keys=\"NAME\")\n", - "region: deepicedrain.Region = deepicedrain.Region.from_gdf(\n", - " gdf=drainage_basins.loc[placename], name=\"Recovery Basin\"\n", - ")\n", + "placename: str = \"siple_coast\" # \"Recovery\" # \"Whillans\"\n", + "try:\n", + " drainage_basins: gpd.GeoDataFrame = drainage_basins.set_index(keys=\"NAME\")\n", + " region: deepicedrain.Region = deepicedrain.Region.from_gdf(\n", + " gdf=drainage_basins.loc[placename], name=\"Recovery Basin\"\n", + " )\n", + "except KeyError:\n", + " pass\n", "df_dhdt: cudf.DataFrame = cudf.read_parquet(\n", " f\"ATLXI/df_dhdt_{placename.lower()}.parquet\"\n", ")" @@ -935,7 +938,7 @@ "source": [ "# Antarctic subglacial lake polygons with EPSG:3031 coordinates\n", "antarctic_lakes: gpd.GeoDataFrame = gpd.read_file(\n", - " filename=\"antarctic_subglacial_lakes.geojson\"\n", + " filename=\"antarctic_subglacial_lakes_3031.geojson\"\n", ")" ] }, @@ -948,17 +951,28 @@ "outputs": [], "source": [ "# Choose one draining/filling lake\n", - "draining: bool = False # False\n", - "placename: str = \"Slessor\" # \"Whillans\"\n", + "draining: bool = False\n", + "placename: str = \"Whillans\" # \"Slessor\" # \"Kamb\" # \"Mercer\" #\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", + "lake = lakes.loc[lakes.inner_dhdt.idxmin() if draining else lakes.inner_dhdt.idxmax()]\n", + "# lake = lakes.query(expr=\"inner_dhdt < 0\" if draining else \"inner_dhdt > 0\").loc[63]\n", "lakedict = {\n", - " 76: \"Subglacial Lake Conway\", # draining lake\n", - " 78: \"Whillans IX\", # filling lake\n", - " 103: \"Slessor 45\", # draining lake\n", - " 108: \"Slessor 23\", # filling lake\n", + " 21: \"Mercer 2b\", # filling lake\n", + " 40: \"Subglacial Lake Conway\", # draining lake\n", + " 48: \"Subglacial Lake Whillans\", # filling lake\n", + " 50: \"Whillans IX\", # filling lake\n", + " 63: \"Kamb 1\", # filling lake\n", + " 65: \"Kamb 12\", # filling lake\n", + " 97: \"MacAyeal 1\", # draining lake\n", + " 109: \"Slessor 45\", # draining lake\n", + " 116: \"Slessor 23\", # filling lake\n", + " 153: \"Recovery IX\", # draining lake\n", + " 157: \"Recovery 3\", # filling lake\n", "}\n", - "region = deepicedrain.Region.from_gdf(gdf=lake, name=lakedict[lake.name])" + "region = deepicedrain.Region.from_gdf(gdf=lake, name=lakedict[lake.name])\n", + "\n", + "print(lake)\n", + "lake.geometry" ] }, { @@ -1064,6 +1078,10 @@ "# 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", " future = client.submit(\n", @@ -1161,14 +1179,20 @@ "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\n", + "# Plot actual track points in green\n", "for track in tracks:\n", - " fig.plot(x=track.x, y=track.y, color=\"green\", style=\"c0.01c\")\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\n", + "# Plot lake boundary in blue\n", "lakex, lakey = lake.geometry.exterior.coords.xy\n", - "fig.plot(x=lakex, y=lakey, pen=\"thin\")\n", + "fig.plot(x=lakex, y=lakey, pen=\"thin,blue,-.\")\n", "# Map frame in kilometre units\n", "fig.basemap(\n", " frame=[\n", @@ -1180,7 +1204,7 @@ " projection=\"X8c\",\n", ")\n", "fig.colorbar(position=\"JMR\", frame=['x+l\"Crossover Error\"', \"y+lm\"])\n", - "fig.savefig(f\"figures/crossover_area_{placename}.png\")\n", + "fig.savefig(f\"figures/{placename}/crossover_area_{placename}_{min_date}_{max_date}.png\")\n", "fig.show()" ] }, @@ -1252,7 +1276,7 @@ "# 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/crossover_point_{placename}_{track1}_{track2}_{min_date}_{max_date}.png\"\n", + " f\"figures/{placename}/crossover_point_{placename}_{track1}_{track2}_{min_date}_{max_date}.png\"\n", ")\n", "fig.show()" ] @@ -1264,8 +1288,8 @@ "outputs": [], "source": [ "# 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 = 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()" ] }, @@ -1280,10 +1304,15 @@ "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", - "fig = deepicedrain.vizplots.plot_crossovers(\n", - " df=df_th, regionname=region.name, elev_var=\"h_norm\"\n", + "fig = deepicedrain.plot_crossovers(\n", + " df=df_th,\n", + " regionname=region.name,\n", + " elev_var=\"h_norm\",\n", + " elev_filter=3 * abs(df.h_X).median(),\n", + ")\n", + "fig.savefig(\n", + " f\"figures/{placename}/crossover_many_normalized_{placename}_{min_date}_{max_date}.png\"\n", ")\n", - "fig.savefig(f\"figures/crossover_many_normalized_{placename}_{min_date}_{max_date}.png\")\n", "fig.show()" ] }, @@ -1315,7 +1344,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.3" + "version": "3.8.6" } }, "nbformat": 4, diff --git a/atlxi_lake.py b/atlxi_lake.py index 29ec9f1..d008d44 100644 --- a/atlxi_lake.py +++ b/atlxi_lake.py @@ -343,11 +343,14 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series: # %% # Save or load dhdt data from Parquet file -placename: str = "Recovery" # "Whillans" -drainage_basins: gpd.GeoDataFrame = drainage_basins.set_index(keys="NAME") -region: deepicedrain.Region = deepicedrain.Region.from_gdf( - gdf=drainage_basins.loc[placename], name="Recovery Basin" -) +placename: str = "siple_coast" # "Recovery" # "Whillans" +try: + drainage_basins: gpd.GeoDataFrame = drainage_basins.set_index(keys="NAME") + region: deepicedrain.Region = deepicedrain.Region.from_gdf( + gdf=drainage_basins.loc[placename], name="Recovery Basin" + ) +except KeyError: + pass df_dhdt: cudf.DataFrame = cudf.read_parquet( f"ATLXI/df_dhdt_{placename.lower()}.parquet" ) @@ -356,23 +359,34 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series: # %% # Antarctic subglacial lake polygons with EPSG:3031 coordinates antarctic_lakes: gpd.GeoDataFrame = gpd.read_file( - filename="antarctic_subglacial_lakes.geojson" + filename="antarctic_subglacial_lakes_3031.geojson" ) # %% # Choose one draining/filling lake -draining: bool = False # False -placename: str = "Slessor" # "Whillans" +draining: bool = False +placename: str = "Whillans" # "Slessor" # "Kamb" # "Mercer" # lakes: gpd.GeoDataFrame = antarctic_lakes.query(expr="basin_name == @placename") -lake = lakes.loc[lakes.maxabsdhdt.idxmin() if draining else lakes.maxabsdhdt.idxmax()] +lake = lakes.loc[lakes.inner_dhdt.idxmin() if draining else lakes.inner_dhdt.idxmax()] +# lake = lakes.query(expr="inner_dhdt < 0" if draining else "inner_dhdt > 0").loc[63] lakedict = { - 76: "Subglacial Lake Conway", # draining lake - 78: "Whillans IX", # filling lake - 103: "Slessor 45", # draining lake - 108: "Slessor 23", # filling lake + 21: "Mercer 2b", # filling lake + 40: "Subglacial Lake Conway", # draining lake + 48: "Subglacial Lake Whillans", # filling lake + 50: "Whillans IX", # filling lake + 63: "Kamb 1", # filling lake + 65: "Kamb 12", # filling lake + 97: "MacAyeal 1", # draining lake + 109: "Slessor 45", # draining lake + 116: "Slessor 23", # filling lake + 153: "Recovery IX", # draining lake + 157: "Recovery 3", # filling lake } region = deepicedrain.Region.from_gdf(gdf=lake, name=lakedict[lake.name]) +print(lake) +lake.geometry + # %% # Subset data to lake of interest placename: str = region.name.lower().replace(" ", "_") @@ -445,6 +459,10 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series: # Parallelized paired crossover analysis futures: list = [] for rgt1, rgt2 in itertools.combinations(rgts, r=2): + # skip if same referencegroundtrack but different laser pair + # as they are parallel and won't cross + if rgt1[:4] == rgt2[:4]: + continue track1 = track_dict[rgt1][["x", "y", "h_corr", "utc_time"]] track2 = track_dict[rgt2][["x", "y", "h_corr", "utc_time"]] future = client.submit( @@ -511,14 +529,20 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series: pygmt.makecpt(cmap="batlow", series=[sumstats[var]["25%"], sumstats[var]["75%"]]) # Map frame in metre units fig.basemap(frame="f", region=plotregion, projection="X8c") -# Plot actual track points +# Plot actual track points in green for track in tracks: - fig.plot(x=track.x, y=track.y, color="green", style="c0.01c") + tracklabel = f"{track.iloc[0].referencegroundtrack} {track.iloc[0].pairtrack}" + fig.plot( + x=track.x, + y=track.y, + pen="thinnest,green,.", + style=f'qN+1:+l"{tracklabel}"+f3p,Helvetica,darkgreen', + ) # Plot crossover point locations fig.plot(x=df.x, y=df.y, color=df.h_X, cmap=True, style="c0.1c", pen="thinnest") -# PLot lake boundary +# Plot lake boundary in blue lakex, lakey = lake.geometry.exterior.coords.xy -fig.plot(x=lakex, y=lakey, pen="thin") +fig.plot(x=lakex, y=lakey, pen="thin,blue,-.") # Map frame in kilometre units fig.basemap( frame=[ @@ -530,7 +554,7 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series: projection="X8c", ) fig.colorbar(position="JMR", frame=['x+l"Crossover Error"', "y+lm"]) -fig.savefig(f"figures/crossover_area_{placename}.png") +fig.savefig(f"figures/{placename}/crossover_area_{placename}_{min_date}_{max_date}.png") fig.show() @@ -587,14 +611,14 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series: # Plot dashed line connecting points fig.plot(x=df_max.t, y=df_max.h, pen=f"faint,blue,-") fig.savefig( - f"figures/crossover_point_{placename}_{track1}_{track2}_{min_date}_{max_date}.png" + f"figures/{placename}/crossover_point_{placename}_{track1}_{track2}_{min_date}_{max_date}.png" ) fig.show() # %% # Plot all crossover height points over time over the lake area -fig = deepicedrain.vizplots.plot_crossovers(df=df_th, regionname=region.name) -fig.savefig(f"figures/crossover_many_{placename}_{min_date}_{max_date}.png") +fig = deepicedrain.plot_crossovers(df=df_th, regionname=region.name) +fig.savefig(f"figures/{placename}/crossover_many_{placename}_{min_date}_{max_date}.png") fig.show() # %% @@ -603,10 +627,15 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series: normfunc = lambda h: h - h.iloc[0] # lambda h: h - h.mean() df_th["h_norm"] = df_th.groupby(by="track1_track2").h.transform(func=normfunc) -fig = deepicedrain.vizplots.plot_crossovers( - df=df_th, regionname=region.name, elev_var="h_norm" +fig = deepicedrain.plot_crossovers( + df=df_th, + regionname=region.name, + elev_var="h_norm", + elev_filter=3 * abs(df.h_X).median(), +) +fig.savefig( + f"figures/{placename}/crossover_many_normalized_{placename}_{min_date}_{max_date}.png" ) -fig.savefig(f"figures/crossover_many_normalized_{placename}_{min_date}_{max_date}.png") fig.show() # %% diff --git a/deepicedrain/vizplots.py b/deepicedrain/vizplots.py index ec29e5c..5458f84 100644 --- a/deepicedrain/vizplots.py +++ b/deepicedrain/vizplots.py @@ -35,7 +35,7 @@ class IceSat2Explorer(param.Parameterized): plot_variable = param.Selector( default="dhdt_slope", objects=["referencegroundtrack", "dhdt_slope", "h_corr"] ) - cycle_number = param.Integer(default=7, bounds=(2, 7)) + cycle_number = param.Integer(default=7, bounds=(2, 8)) dhdt_range = param.Range(default=(1.0, 10.0), bounds=(0.0, 20.0)) rasterize = param.Boolean(default=True) datashade = param.Boolean(default=False) @@ -271,6 +271,7 @@ def plot_crossovers( time_var: str = "t", track_var: str = "track1_track2", spacing: float = 2.5, + elev_filter: float = 1.0, ) -> pygmt.Figure: """ Plot to show how elevation is changing at many crossover points over time. @@ -316,6 +317,9 @@ def plot_crossovers( Provide as a 'dy' increment, this is passed on to `pygmt.info` and used to round up and down the y axis (elev_var) limits for a nicer plot frame. Default is 2.5. + elev_filter : float + Minimum elevation change required for the crossover point to show up + on the plot. Default is 1.0 (metres). Returns ------- @@ -359,7 +363,7 @@ def plot_crossovers( ): df_ = df.loc[indexes].sort_values(by=time_var) # plot only > 1 metre height change - if df_[elev_var].max() - df_[elev_var].min() > 1.0: + if df_[elev_var].max() - df_[elev_var].min() > elev_filter: track1, track2 = track1_track2.split("x") fig.plot( x=df_[time_var], @@ -374,5 +378,6 @@ def plot_crossovers( fig.plot( x=df_[time_var], y=df_[elev_var], Z=i, pen=f"faint,+z,-", cmap=True ) - fig.legend(S=0.5, position="JMR+JMR+o0.2c", box="+gwhite+p1p") + with pygmt.config(FONT_ANNOT_PRIMARY="9p"): + fig.legend(S=0.8, position="JTR+jTL+o0.2c", box="+gwhite+p1p") return fig From 38237d10b9e092b916eb1578a119b15eedbb45da Mon Sep 17 00:00:00 2001 From: Wei Ji Date: Mon, 26 Oct 2020 22:57:20 +1300 Subject: [PATCH 5/5] :dizzy: 3D view of interpolated ice surface elev over ICESat-2 cycles Visualizing ice surface changes in 3D, because it's more intuitive to the eye! There's a bit of an art to 3D perspective plots, and PyGMT's `grdview` does a fine job at it, though there's many lines of code involved! The figure consists of two plots. On the left is a shaded ice surface elevation grid overlaid with the subglacial lake outline in yellow. On the right is an elevation difference/anomaly grid view (relative to ICESat-2 Cycle 3). Still need to wait for `subplot` and `plot3d` to be wrapped in PyGMT to simplify the code and make things look more polished. Animated GIF is made using imagemagick's `convert` in bash as it was higher quality (and less lines of code) than Pillow. Also have a code cell showing an alternative HvPlot-based 2D map of the same time-series grid, complete with an interactive slider to cycle through each ICESat-2 cycle. --- atlxi_lake.ipynb | 216 +++++++++++++++++++++++++++++++++++++++++++++-- atlxi_lake.py | 122 +++++++++++++++++++++++++- 2 files changed, 331 insertions(+), 7 deletions(-) diff --git a/atlxi_lake.ipynb b/atlxi_lake.ipynb index 2e8bdb6..f4aa4d7 100644 --- a/atlxi_lake.ipynb +++ b/atlxi_lake.ipynb @@ -42,6 +42,7 @@ "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", @@ -49,6 +50,7 @@ "import scipy.spatial\n", "import shapely.geometry\n", "import tqdm\n", + "import xarray as xr\n", "import zarr\n", "\n", "import deepicedrain" @@ -905,12 +907,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Select a lake to examine" + "# Select a subglacial lake to examine" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": { "lines_to_next_cell": 2 }, @@ -932,7 +934,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ @@ -944,11 +946,41 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "metadata": { "lines_to_next_cell": 1 }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "basin_name Whillans\n", + "num_points 3276\n", + "outer_dhdt 0.377611\n", + "outer_std 1.42712\n", + "outer_mad 0.083278\n", + "inner_dhdt 2.12607\n", + "maxabsdhdt 8.65662\n", + "refgtracks 74|135|196|266|327|388|516|577|638|769|830|101...\n", + "geometry POLYGON ((-444681.0351994165 -545210.454471163...\n", + "Name: 50, dtype: object\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "# Choose one draining/filling lake\n", "draining: bool = False\n", @@ -977,7 +1009,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, "metadata": { "lines_to_next_cell": 2 }, @@ -988,6 +1020,178 @@ "df_lake: cudf.DataFrame = region.subset(data=df_dhdt)" ] }, + { + "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)\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", + "grid_region = pygmt.info(table=df_lake[[\"x\", \"y\"]].to_pandas(), spacing=\"s250\")\n", + "z_limits: tuple = (float(ds_lake.z.min()), float(ds_lake.z.max())) # original z limits\n", + "grid_region: np.ndarray = np.append(arr=grid_region, values=z_limits)\n", + "\n", + "ds_lake_norm: xr.Dataset = ds_lake - ds_lake.sel(cycle_number=3).z\n", + "z_norm_limits: tuple = (float(ds_lake_norm.z.min()), float(ds_lake_norm.z.max()))\n", + "grid_region_norm: np.ndarray = np.append(arr=grid_region[:4], values=z_norm_limits)\n", + "\n", + "print(f\"Elevation limits are: {z_limits}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 3D plots of gridded ice surface elevation over time\n", + "azimuth: float = 157.5 # 202.5 # 270\n", + "elevation: float = 45 # 60\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", + " grdview_kwargs = dict(\n", + " cmap=True,\n", + " zscale=0.1, # zscaling factor, default to 10x vertical exaggeration\n", + " surftype=\"sim\", # surface, image and mesh plot\n", + " perspective=[azimuth, elevation], # perspective using azimuth/elevation\n", + " # W=\"c0.05p,black,solid\", # draw contours\n", + " )\n", + "\n", + " fig = pygmt.Figure()\n", + " # grid = ds_lake.sel(cycle_number=cycle).z\n", + " grid = f\"figures/{placename}/h_corr_{placename}_cycle_{cycle}.nc\"\n", + " pygmt.makecpt(cmap=\"lapaz\", series=z_limits)\n", + " fig.grdview(\n", + " grid=grid,\n", + " projection=\"X10c\",\n", + " region=grid_region,\n", + " shading=True,\n", + " frame=[\n", + " f'SWZ+t\"{region.name}\"', # plot South, West axes, and Z-axis\n", + " 'xaf+l\"Polar Stereographic X (m)\"', # add x-axis annotations and minor ticks\n", + " 'yaf+l\"Polar Stereographic Y (m)\"', # add y-axis annotations and minor ticks\n", + " f'zaf+l\"Elevation (m)\"', # add z-axis annotations, minor ticks and axis label\n", + " ],\n", + " **grdview_kwargs,\n", + " )\n", + "\n", + " # Plot lake boundary outline\n", + " # TODO wait for plot3d to plot lake boundary points at correct height\n", + " df = pd.DataFrame([region.bounds()]).values\n", + " points = pd.DataFrame(\n", + " data=[point for point in lake.geometry.exterior.coords], columns=(\"x\", \"y\")\n", + " )\n", + " df_xyz = pygmt.grdtrack(points=points, grid=grid, newcolname=\"z\")\n", + " fig.plot(\n", + " data=df_xyz.values,\n", + " region=grid_region,\n", + " pen=\"1.5p,yellow2\",\n", + " Jz=True, # zscale\n", + " p=f\"{azimuth}/{elevation}/{df_xyz.z.median()}\", # perspective\n", + " # label='\"Subglacial Lake X\"'\n", + " )\n", + "\n", + " # Plot normalized elevation change\n", + " grid = ds_lake_norm.sel(cycle_number=cycle).z\n", + " if cycle == 3:\n", + " # add some tiny random noise to make plot work\n", + " grid = grid + np.random.normal(scale=1e-8, size=grid.shape)\n", + " pygmt.makecpt(cmap=\"vik\", series=z_norm_limits)\n", + " fig.grdview(\n", + " grid=grid,\n", + " region=grid_region_norm,\n", + " frame=[\n", + " f'SEZ2+t\"Cycle {cycle} at {time_sec}\"', # plot South, East axes, and Z-axis\n", + " 'xaf+l\"Polar Stereographic X (m)\"', # add x-axis annotations and minor ticks\n", + " 'yaf+l\"Polar Stereographic Y (m)\"', # add y-axis annotations and minor ticks\n", + " f'zaf+l\"Elev Change (m)\"', # add z-axis annotations, minor ticks and axis label\n", + " ],\n", + " X=\"10c\", # xshift\n", + " **grdview_kwargs,\n", + " )\n", + "\n", + " fig.savefig(f\"figures/{placename}/dsm_{placename}_cycle_{cycle}.png\")\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Make a animated GIF of changing ice surface from the PNG files\n", + "gif_fname: str = (\n", + " f\"figures/{placename}/dsm_{placename}_cycles_{cycles[0]}-{cycles[-1]}.gif\"\n", + ")\n", + "# !convert -delay 120 -loop 0 figures/{placename}/dsm_*.png {gif_fname}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "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, diff --git a/atlxi_lake.py b/atlxi_lake.py index d008d44..edebad3 100644 --- a/atlxi_lake.py +++ b/atlxi_lake.py @@ -44,6 +44,7 @@ import dask import dask.array import geopandas as gpd +import hvplot.xarray import numpy as np import pandas as pd import panel as pn @@ -51,6 +52,7 @@ import scipy.spatial import shapely.geometry import tqdm +import xarray as xr import zarr import deepicedrain @@ -339,7 +341,7 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series: # %% [markdown] -# # Select a lake to examine +# # Select a subglacial lake to examine # %% # Save or load dhdt data from Parquet file @@ -393,6 +395,124 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series: df_lake: cudf.DataFrame = region.subset(data=df_dhdt) +# %% [markdown] +# ## Create an interpolated ice surface elevation grid for each ICESat-2 cycle + +# %% +# Generate gridded time-series of ice elevation over lake +cycles: tuple = (3, 4, 5, 6, 7, 8) +os.makedirs(name=f"figures/{placename}", exist_ok=True) +ds_lake: xr.Dataset = deepicedrain.spatiotemporal_cube( + table=df_lake.to_pandas(), + placename=placename, + cycles=cycles, + folder=f"figures/{placename}", +) +ds_lake.to_netcdf(path=f"figures/{placename}/xyht_{placename}.nc", mode="w") + +# %% +# Get 3D grid_region (xmin/xmax/ymin/ymax/zmin/zmax), +# and calculate normalized z-values as Elevation delta relative to Cycle 3 +grid_region = pygmt.info(table=df_lake[["x", "y"]].to_pandas(), spacing="s250") +z_limits: tuple = (float(ds_lake.z.min()), float(ds_lake.z.max())) # original z limits +grid_region: np.ndarray = np.append(arr=grid_region, values=z_limits) + +ds_lake_norm: xr.Dataset = ds_lake - ds_lake.sel(cycle_number=3).z +z_norm_limits: tuple = (float(ds_lake_norm.z.min()), float(ds_lake_norm.z.max())) +grid_region_norm: np.ndarray = np.append(arr=grid_region[:4], values=z_norm_limits) + +print(f"Elevation limits are: {z_limits}") + +# %% +# 3D plots of gridded ice surface elevation over time +azimuth: float = 157.5 # 202.5 # 270 +elevation: float = 45 # 60 +for cycle in tqdm.tqdm(iterable=cycles): + time_nsec: pd.Timestamp = df_lake[f"utc_time_{cycle}"].to_pandas().mean() + time_sec: str = np.datetime_as_string(arr=time_nsec.to_datetime64(), unit="s") + + grdview_kwargs = dict( + cmap=True, + zscale=0.1, # zscaling factor, default to 10x vertical exaggeration + surftype="sim", # surface, image and mesh plot + perspective=[azimuth, elevation], # perspective using azimuth/elevation + # W="c0.05p,black,solid", # draw contours + ) + + fig = pygmt.Figure() + # grid = ds_lake.sel(cycle_number=cycle).z + grid = f"figures/{placename}/h_corr_{placename}_cycle_{cycle}.nc" + pygmt.makecpt(cmap="lapaz", series=z_limits) + fig.grdview( + grid=grid, + projection="X10c", + region=grid_region, + shading=True, + frame=[ + f'SWZ+t"{region.name}"', # plot South, West axes, and Z-axis + 'xaf+l"Polar Stereographic X (m)"', # add x-axis annotations and minor ticks + 'yaf+l"Polar Stereographic Y (m)"', # add y-axis annotations and minor ticks + f'zaf+l"Elevation (m)"', # add z-axis annotations, minor ticks and axis label + ], + **grdview_kwargs, + ) + + # Plot lake boundary outline + # TODO wait for plot3d to plot lake boundary points at correct height + df = pd.DataFrame([region.bounds()]).values + points = pd.DataFrame( + data=[point for point in lake.geometry.exterior.coords], columns=("x", "y") + ) + df_xyz = pygmt.grdtrack(points=points, grid=grid, newcolname="z") + fig.plot( + data=df_xyz.values, + region=grid_region, + pen="1.5p,yellow2", + Jz=True, # zscale + p=f"{azimuth}/{elevation}/{df_xyz.z.median()}", # perspective + # label='"Subglacial Lake X"' + ) + + # Plot normalized elevation change + grid = ds_lake_norm.sel(cycle_number=cycle).z + if cycle == 3: + # add some tiny random noise to make plot work + grid = grid + np.random.normal(scale=1e-8, size=grid.shape) + pygmt.makecpt(cmap="vik", series=z_norm_limits) + fig.grdview( + grid=grid, + region=grid_region_norm, + frame=[ + f'SEZ2+t"Cycle {cycle} at {time_sec}"', # plot South, East axes, and Z-axis + 'xaf+l"Polar Stereographic X (m)"', # add x-axis annotations and minor ticks + 'yaf+l"Polar Stereographic Y (m)"', # add y-axis annotations and minor ticks + f'zaf+l"Elev Change (m)"', # add z-axis annotations, minor ticks and axis label + ], + X="10c", # xshift + **grdview_kwargs, + ) + + fig.savefig(f"figures/{placename}/dsm_{placename}_cycle_{cycle}.png") +fig.show() + +# %% +# Make a animated GIF of changing ice surface from the PNG files +gif_fname: str = ( + f"figures/{placename}/dsm_{placename}_cycles_{cycles[0]}-{cycles[-1]}.gif" +) +# !convert -delay 120 -loop 0 figures/{placename}/dsm_*.png {gif_fname} + +# %% +# HvPlot 2D interactive view of ice surface elevation grids over each ICESat-2 cycle +dashboard: pn.layout.Column = pn.Column( + ds_lake.hvplot.image(x="x", y="y", clim=z_limits, cmap="gist_earth", data_aspect=1) + # * ds_lake.hvplot.contour(x="x", y="y", clim=z_limits, data_aspect=1) +) +dashboard.show(port=30227) + +# %% [markdown] +# ## Along track plots of ice surface elevation change over time + # %% # Select a few Reference Ground tracks to look at rgts: list = [int(rgt) for rgt in lake.refgtracks.split("|")]