KangerFjord.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "89d33882-25fb-4204-9d1b-bbc0a3ba5c1b",
   "metadata": {},
   "source": [
    "# Demonstrate `eulerian.Position` object with Fjord\n",
    "\n",
    "Author: Wenrui Jiang, 14 June 2023\n",
    "\n",
    "The `eulerian.Position` object is really what connects a point and the numerical model. Its `interpolate` method really is the core of this package. We're going to use a rather interesting example to demonstrate the functionalities of `eulerian.Position`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8fbebdd-e8c7-4308-8c38-06d900bbd216",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import oceanspy as ospy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import seaduck as sd\n",
    "from functools import partial\n",
    "import xarray as xr\n",
    "import cmocean\n",
    "\n",
    "mpl.rcParams[\"figure.dpi\"] = 300"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85a2e829-711c-49e4-8af8-b758f2942d3d",
   "metadata": {},
   "source": [
    "Here on [SciServer](https://sciserver.org/https://sciserver.org/), we have an interesting dataset simulating the interaction between background circulation and Kangerdlugssuaq Fjord. More information can be found below, and see the paper by [Fraser et al., 2018](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JC014435https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JC014435):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63546d43-f4e7-4817-9217-511de64da083",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fjord = ospy.open_oceandataset.from_catalog(\"KangerFjord\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ac50b0d-1fea-4c3d-8c6e-eae3aa654479",
   "metadata": {},
   "source": [
    "```\n",
    "Opening KangerFjord.\n",
    "A realistic numerical model constructed to simulate the oceanic conditions\n",
    "and circulation in a large southeast Greenland fjord (Kangerdlugssuaq) and\n",
    "the adjacent shelf sea region during winter 2007–2008.\n",
    "Citation:\n",
    " * Fraser et al., 2018 - JGR.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43d2fea7-a78e-4e99-be4a-cdcb113009cc",
   "metadata": {},
   "source": [
    "Let's first explore the dataset a bit. We are going to use [ETOPO](https://www.ncei.noaa.gov/products/etopo-global-relief-modelhttps://www.ncei.noaa.gov/products/etopo-global-relief-model) dataset to give you an idea where this domain is located."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2192cf7a-e32e-4368-9885-3109666bc468",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "etopo = ospy.open_oceandataset.from_catalog(\"ETOPO\")\n",
    "etopo = etopo._ds.sel(X=slice(-40.6, -13.1), Y=slice(62.4, 70.5)).where(\n",
    "    etopo._ds.Depth > 0\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "474e3c37-8776-4c5e-82ca-942e1344330b",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "vmax = 0.2\n",
    "cmap = cmocean.cm.balance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e04dba6e-7b17-4766-8e8d-5a88e7cff160",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "c = plt.pcolormesh(\n",
    "    fjord._ds.XC, fjord._ds.YC, fjord._ds.Eta[263], vmax=vmax, vmin=-vmax, cmap=cmap\n",
    ")\n",
    "plt.pcolormesh(\n",
    "    etopo.X,\n",
    "    etopo.Y,\n",
    "    etopo.Depth,\n",
    "    cmap=\"bone\",\n",
    "    vmin=-3000,\n",
    "    vmax=4500,\n",
    "    zorder=20,\n",
    ")\n",
    "plt.gca().set_facecolor(\"lightsteelblue\")\n",
    "plt.title(\"Model domain and the sea surface height\")\n",
    "plt.xlabel(\"Longitude\")\n",
    "plt.ylabel(\"Latitude\")\n",
    "plt.colorbar(c, label=\"m\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c82117b5-159f-4eb1-9282-d125851b4bdf",
   "metadata": {
    "tags": [
     "mdformat-skip"
    ]
   },
   "source": [
    "## Doing Interpolation\n",
    "We are going to use the sea surface height field $\\eta$ as an example. \n",
    "\n",
    "First, we are going to convert the `xarray.Dataset` to `seaduck.OceData`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c77051ef-121d-453c-8066-402629e58598",
   "metadata": {},
   "outputs": [],
   "source": [
    "tub = sd.OceData(fjord._ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58e0bf9a-bc27-432f-9b90-1e24892ac470",
   "metadata": {},
   "source": [
    "For now, the interpolation points are simply the gridpoints of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "983b4b44-0120-407b-93bb-e03c8a24aa8a",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "x = tub.XC.ravel()\n",
    "y = tub.YC.ravel()\n",
    "z = np.ones_like(x) * -5.0\n",
    "t = np.ones_like(x) * sd.utils.convert_time(\"2008-01-01\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51ba7082-52e2-418d-9be6-8db68883271f",
   "metadata": {},
   "source": [
    "Let's create the `eulerian.Position` object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88a41587-2ef8-4ee8-8d91-23f8c8b1381c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "p = sd.Position().from_latlon(x, y, z, t, data=tub)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2567e4d-52d5-4c99-b5e7-7c36fa08d1ac",
   "metadata": {},
   "source": [
    "Two interpolation kernels are defined here, both of them are default kernels used by the package. However, we are going to \"hack\" one of them to demonstrate the \"cascade\" capacity of interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6ae57d1-174d-43b0-a423-4a688a0a0a71",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "kernel = sd.KnW()\n",
    "kernel_to_be_hacked = sd.KnW()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2c8b35e-cfb0-43e7-80f7-2f2de89d1cbf",
   "metadata": {
    "tags": [
     "mdformat-skip"
    ]
   },
   "source": [
    "First, we do the normal interpolation on $\\eta$ in the normal way, and plot it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f9e7fdb-d8a0-4b8e-a002-ca639e10b154",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "eta = p.interpolate(\"Eta\", kernel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b1ca767-e44c-4acd-a9f9-cc16bbf49270",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "plt.scatter(x, y, c=eta, s=0.7, vmax=vmax, vmin=-vmax, cmap=cmap)\n",
    "plt.title(\"Interpolated sea surface height\")\n",
    "plt.colorbar(label=\"m\")\n",
    "plt.xlabel(\"Longitude\")\n",
    "plt.ylabel(\"Latitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd4d9629-eb46-4f2e-a54b-31e07fa58176",
   "metadata": {},
   "source": [
    "Now let's \"hack\" the code to make it return the size of the kernels used.\n",
    "\n",
    "The details of the hack is not very important. But if you're interested you can read the inline comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3051c0b-df9f-49dc-8fbd-e5a48d01c4fc",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "# The partials of this function is going to replace all the weight functions.\n",
    "# For legal weight functions, the sum of contribution (weight) should be one,\n",
    "# but here every neighboring point has contribution of one\n",
    "\n",
    "\n",
    "def dummy_weight(rx, ry, ans):\n",
    "    n = len(rx)\n",
    "    return np.ones((n, ans))\n",
    "\n",
    "\n",
    "# This is the step that overwrites the existing weight function. Highly not recommended.\n",
    "kernel_to_be_hacked.hfuncs = [partial(dummy_weight, ans=len(k)) for k in kernel.kernels]\n",
    "# Create a dataset of ones with the same shape as Eta.\n",
    "tub[\"ones\"] = xr.ones_like(tub[\"Eta\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "824b9580-1f54-4476-8c0f-d0db38c1a34a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "how_many = p.interpolate(\"ones\", kernel_to_be_hacked)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "987d48f7-1f8a-4ca1-965a-553050227c3e",
   "metadata": {},
   "source": [
    "Now, we can look at what the interpolation package does under the hood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86e135f9-364a-4176-8902-bce7f3770b43",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "plt.scatter(x, y, c=how_many, s=0.7, cmap=\"Set1\", vmax=9.5, vmin=0.5)\n",
    "plt.title(\"The number of grid points used for interpolation\")\n",
    "plt.colorbar(label=\"Point count\")\n",
    "plt.xlabel(\"Longitude\")\n",
    "plt.ylabel(\"Latitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9007166-2d6f-4d78-bc1f-c602fcafa371",
   "metadata": {},
   "source": [
    "## Filling between\n",
    "\n",
    "As you have seen earlier, the grid of this dataset has very uneven spacing.\n",
    "\n",
    "Well, the strength of the `eulerian.Position` object is filling in information (interpolation) between the grid points.\n",
    "So let's do that. We specify part of the model domain to fill in using the parameters to the `utils.easy_3d_cube` call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42b1d75d-0751-4f6c-a6fe-f4a88d5b3b12",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "Nlon = 600\n",
    "Nlat = 100\n",
    "Ndep = 1\n",
    "ax, ay, az, at = sd.utils.easy_3d_cube(\n",
    "    (-34.5, -28.5, Nlon), (66.5, 67.0, Nlat), (-5.0, -5.0, Ndep), \"2008-01-01\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7b55c45-2528-4c0e-bbbc-51a8906f6e1b",
   "metadata": {},
   "source": [
    "Create `eulerian.Position` in between, and plot it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "234e99c9-6c80-4190-a262-4fde0ffdf737",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "addition = sd.Position().from_latlon(ax, ay, az, at, data=tub)\n",
    "more_eta = addition.interpolate(\"Eta\", kernel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "035b10f6-b3fc-41f8-aaef-02c8fce51eb4",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "plt.scatter(x, y, c=eta, s=0.7, vmax=vmax, vmin=-vmax, cmap=cmap)\n",
    "plt.scatter(ax, ay, c=more_eta, s=0.7, vmax=vmax, vmin=-vmax, cmap=cmap)\n",
    "plt.title(\"Interpolated sea surface height (Reinforced)\")\n",
    "plt.xlabel(\"Longitude\")\n",
    "plt.ylabel(\"Latitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "541b25d2-388c-4c1e-8e71-1fb6bd9d091a",
   "metadata": {},
   "source": [
    "I'd say the filling-in is done pretty well!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Oceanography",
   "language": "python",
   "name": "oceanography"
  },
  "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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}