Notebook to implement DEM smoothing.

notebooks/dem_smoothing.ipynb

@@ -0,0 +1,353 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "5b7dd5df-edd1-4e40-b83a-7ee663146c29",
+   "metadata": {},
+   "source": [
+    "# Ice-thickness length scale smoothing\n",
+    "\n",
+    "We are implmenting eq 2 and 3 from https://polarresearch.net/index.php/polar/article/view/3498/9172"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b1c1bebf-42cd-48cf-a7bb-669237385cf2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.ndimage import gaussian_filter, generic_filter\n",
+    "import xarray as xr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b9dcb1d1-aa7a-4308-ad96-7b706c7ffc8e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pathlib import Path\n",
+    "from typing import Dict, List, Optional, Union\n",
+    "\n",
+    "import cartopy.crs as ccrs\n",
+    "import cf_xarray.units  # pylint: disable=unused-import\n",
+    "import geopandas as gp\n",
+    "import matplotlib\n",
+    "import numpy as np\n",
+    "import pint_xarray  # pylint: disable=unused-import\n",
+    "import pylab as plt\n",
+    "import seaborn as sns\n",
+    "import xarray as xr\n",
+    "from matplotlib import cm, colors\n",
+    "from matplotlib.colors import LightSource\n",
+    "from shapely import get_coordinates\n",
+    "\n",
+    "from glacier_flow_tools.utils import (\n",
+    "    blend_multiply,\n",
+    "    figure_extent,\n",
+    "    get_dataarray_extent,\n",
+    "    register_colormaps,\n",
+    ")\n",
+    "\n",
+    "register_colormaps()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c2d14042-3ef9-404f-81b6-4ce0af8bc28a",
+   "metadata": {},
+   "source": [
+    "## This plotting function makes a nice hillshade so it's easier to validate the results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c33add49-a58a-45b4-85e0-17f86ec8a200",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_glacier(\n",
+    "    surface: xr.DataArray,\n",
+    "    overlay: xr.DataArray,\n",
+    "    cmap: str = \"viridis\",\n",
+    "    vmin: float = 10,\n",
+    "    vmax: float = 1500,\n",
+    "    ticks: Union[List[float], np.ndarray] = [10, 100, 250, 500, 750, 1500],\n",
+    "    fontsize: float = 6,\n",
+    "    figwidth: float = 3.2,\n",
+    "):\n",
+    "    \"\"\"\n",
+    "    Plot a surface over a hillshade, add profile and correlation coefficient.\n",
+    "\n",
+    "    This function plots a surface over a hillshade, adds a profile and correlation coefficient.\n",
+    "    The plot is saved as a PDF file in the specified result directory.\n",
+    "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    surface : xr.DataArray\n",
+    "        The surface to be plotted over the hillshade.\n",
+    "    overlay : xr.DataArray\n",
+    "        The overlay to be added to the plot.\n",
+    "    result_dir : Union[str, Path]\n",
+    "        The directory where the result PDF file will be saved.\n",
+    "    cmap : str, optional\n",
+    "        The colormap to be used for the plot, by default \"viridis\".\n",
+    "    vmin : float, optional\n",
+    "        The minimum value for the colormap, by default 10.\n",
+    "    vmax : float, optional\n",
+    "        The maximum value for the colormap, by default 1500.\n",
+    "    ticks : Union[List[float], np.ndarray], optional\n",
+    "        The ticks to be used for the colorbar, by default [10, 100, 250, 500, 750, 1500].\n",
+    "    fontsize : float, optional\n",
+    "        The font size to be used for the plot, by default 6.\n",
+    "    figwidth : float, optional\n",
+    "        The width of the figure in inches, by default 3.2.\n",
+    "\n",
+    "    Examples\n",
+    "    --------\n",
+    "    >>> plot_glacier(profile_series, surface, overlay, '/path/to/result_dir')\n",
+    "    \"\"\"\n",
+    "    plt.rcParams[\"font.size\"] = fontsize\n",
+    "    cartopy_crs = ccrs.NorthPolarStereo(central_longitude=-45, true_scale_latitude=70, globe=None)\n",
+    "    # Shade from the northwest, with the sun 45 degrees from horizontal\n",
+    "    light_source = LightSource(azdeg=315, altdeg=45)\n",
+    "    glacier_overlay = overlay\n",
+    "    glacier_surface = surface.interp_like(glacier_overlay)\n",
+    "\n",
+    "    extent = get_dataarray_extent(glacier_overlay)\n",
+    "    norm = colors.Normalize(vmin=vmin, vmax=vmax)\n",
+    "    mapper = cm.ScalarMappable(norm=norm, cmap=cmap)\n",
+    "\n",
+    "    v = mapper.to_rgba(glacier_overlay.to_numpy())\n",
+    "    z = glacier_surface.to_numpy()\n",
+    "\n",
+    "    ar = 1.0  # initial aspect ratio for first trial\n",
+    "    wi = figwidth  # width in inches\n",
+    "    hi = wi * ar  # height in inches\n",
+    "\n",
+    "    fig = plt.figure(figsize=(wi, hi))\n",
+    "    ax = fig.add_subplot(111, projection=cartopy_crs)\n",
+    "    rgb = light_source.shade_rgb(v, elevation=z, vert_exag=0.01, blend_mode=blend_multiply)\n",
+    "    # Use a proxy artist for the colorbar...\n",
+    "    im = ax.imshow(v, cmap=cmap, vmin=vmin, vmax=vmax)\n",
+    "    im.remove()\n",
+    "    ax.imshow(rgb, extent=extent, origin=\"upper\", transform=cartopy_crs)\n",
+    "    ax.gridlines(\n",
+    "        draw_labels={\"top\": \"x\", \"left\": \"y\"},\n",
+    "        dms=True,\n",
+    "        xlocs=np.arange(-50, 0, 1),\n",
+    "        ylocs=np.arange(50, 88, 1),\n",
+    "        x_inline=False,\n",
+    "        y_inline=False,\n",
+    "        rotate_labels=20,\n",
+    "        ls=\"dotted\",\n",
+    "        color=\"k\",\n",
+    "    )\n",
+    "\n",
+    "    fig.colorbar(im, ax=ax, shrink=0.5, pad=0.025, label=overlay.units, extend=\"both\", ticks=ticks)\n",
+    "    plt.draw()\n",
+    "\n",
+    "    # Get proper ratio here\n",
+    "    xmin, xmax = ax.get_xbound()\n",
+    "    ymin, ymax = ax.get_ybound()\n",
+    "    y2x_ratio = (ymax - ymin) / (xmax - xmin)\n",
+    "    fig.set_figheight(wi * y2x_ratio)\n",
+    "    fig.tight_layout()\n",
+    "    plt.close()\n",
+    "    return fig\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "69a1e251-5cd8-438f-bc7b-954be990aae1",
+   "metadata": {},
+   "source": [
+    "## We use Bedmachine as it has both surface and ice thickness"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4898c5ff-0226-4e6f-8e9a-abf7bccba346",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = xr.open_dataset(\"/Users/andy/Google Drive/My Drive/data/MCdataset/BedMachineGreenland-v5.nc\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b7d2aadc-b822-4fe5-bbc6-d476ecf8a469",
+   "metadata": {},
+   "source": [
+    "## Select Jakobshavn"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9b030608-789e-4a3c-b189-2846fe95ea72",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "jak_ds = ds.sel(x=slice(-226_000, -140_000), y=slice(-2_250_000, -2_300_000))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0fde9da1-50db-4a6a-ad59-4ef9521d9c14",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dem = jak_ds[\"surface\"]\n",
+    "plot_glacier(dem, dem, cmap=\"Grays\", figwidth=16)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e1ecd8b3-b65c-4f9f-a236-a7763a45d583",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dem_smoothed_gaussian = jak_ds[\"surface\"].copy()\n",
+    "dem_smoothed_gaussian.values = gaussian_filter(jak_ds[\"surface\"], 2)\n",
+    "plot_glacier(dem_smoothed_gaussian, dem_smoothed_gaussian, cmap=\"Grays\", figwidth=16)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a0cc72c2-e58f-4575-b6de-18548b2a96db",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "H = jak_ds[\"thickness\"].copy()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8c595f78-0912-4df7-b152-fd637eb2aa73",
+   "metadata": {},
+   "source": [
+    "## Brinkerhoff, Aschwanden and Fahnestock (2021): DOI:10.1017/jog.2020.112\n",
+    "\n",
+    "I think Eq 45 should give us a distance scaled sigma?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3a076182-0572-4c7c-b18e-1a93f3ad1295",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "l = 4  # length scale\n",
+    "L = l * (H.to_numpy() @ H.to_numpy().T)**(1/2) # I am not sure this is correct"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bbad2949-f9fd-4fcc-8f1d-cfa09db65169",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.spatial.distance import pdist, squareform"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "662abbdc-8788-464b-82c7-1d4807231baf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "X, Y = np.meshgrid(H.x, H.y)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "16730aba-8771-40e5-8c74-a9c9ad01cea4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "D = squareform(pdist(Y))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "df3e8c71-3fdb-4023-8ad4-127dcfd9cd6e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Sigma = (1 + D**2 / (2 * L**2))**(-1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bdac6704-0089-4413-b26a-a003699786a1",
+   "metadata": {},
+   "source": [
+    "## Gaussian smoothing\n",
+    "\n",
+    "$$ \\omega=\\frac{1}{2 \\pi \\sigma^2} e^{\\frac{x^2+y^2}{2 \\sigma^2}}$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9fe11450-6867-4d1a-9577-edb1a8210cab",
+   "metadata": {},
+   "source": [
+    "How do we implment this as a **generic_filter**? The implementation of **gaussian_filter** is here: https://github.com/scipy/scipy/blob/v1.13.1/scipy/ndimage/_filters.py#L286-L390.\n",
+    "It uses *np.correlate* instead of convolve, and loops over all axis. Can we do that too?\n",
+    "\n",
+    "Maybe we can use https://github.com/scipy/scipy/blob/v1.13.1/scipy/ndimage/_filters.py#L186 for inspiration. We need to compute $\\sigma$ scaled by ice thickness in a sensible way. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2dca7e41-2e62-4576-8996-13b25bafa88d",
+   "metadata": {},
+   "source": [
+    "## Triangular smoothing\n",
+    "\n",
+    "$$ \\omega = \\max\\left( \\frac{\\sqrt{x^2+y^2}}{\\sigma},0 \\right)$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1af63265-8c5b-4e7e-b4fa-86cf6803d122",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}