ECCO_plot_stations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "300645b1",
   "metadata": {},
   "source": [
    "# `OceInterp` interpolation in the ECCO global state estimate.\n",
    "\n",
    "Extract temperature/salinity profiles at specified longitudes, latitudes, and times. This process mimics hydrographic stations from in-situ measurements (like a ship CTD station, or Argo float profile).\n",
    "\n",
    "This notebook uses [Oceanspy](https://oceanspy.readthedocs.io/en/latest/) and demonstrates the interface to the Poseidon-viewer on SciServer.\n",
    "\n",
    "Author: Tom Haine & Wenrui Jiang, Jun '23"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61c4c3ec",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import seaduck as sd\n",
    "import oceanspy as ospy\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = 12, 8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b048a92",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-05T20:00:52.452205Z",
     "iopub.status.busy": "2023-07-05T20:00:52.451627Z",
     "iopub.status.idle": "2023-07-05T20:00:52.481674Z",
     "shell.execute_reply": "2023-07-05T20:00:52.476797Z",
     "shell.execute_reply.started": "2023-07-05T20:00:52.452149Z"
    }
   },
   "source": [
    "```{admonition} Access ECCO\n",
    "The global MITgcm run is the ECCO state estimate [(Forget et al, 2015)](https://gmd.copernicus.org/articles/8/3071/2015/) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. \n",
    "\n",
    "`ecco = ospy.open_oceandataset.from_catalog(\"ECCO\")`\n",
    "\n",
    "Click [here](https://dev-poseidon-ocean.pantheonsite.io/products/datasets/) for a full list of the dataset hosted and [here](https://oceanspy.readthedocs.io/en/latest/datasets.html#igpwinterhttps://oceanspy.readthedocs.io/en/latest/datasets.html#igpwinter) to find out more.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b78b5ef8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "od = ospy.open_oceandataset.from_catalog(\"ECCO\")\n",
    "\n",
    "od._ds = od._ds.drop_vars({\"k\", \"k_u\", \"k_p1\", \"k_l\"})\n",
    "co_list = [var for var in od._ds.variables if \"time\" not in od._ds[var].dims]\n",
    "od._ds = od._ds.set_coords(co_list)\n",
    "\n",
    "od._ds = od._ds.drop_vars(\"time_bnds\")\n",
    "od._ds[\"Temp\"] = od._ds[\"THETA\"]\n",
    "od._ds[\"S\"] = od._ds[\"SALT\"]\n",
    "this_time = \"1992-03-16T12:00:00.000000000\"\n",
    "# Select time of interest.\n",
    "# This is cut and pasted (with editing to convert format) from Poseidon-viewer,\n",
    "# but in future should be returned in the JSON object.\n",
    "\n",
    "varList = [\"Temp\", \"S\"]  # Select variables of interest.\n",
    "Nvars = len(varList)\n",
    "od._ds = od._ds.drop_vars([var for var in od._ds.data_vars if var not in varList])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6314e77e",
   "metadata": {},
   "source": [
    "Define coordinates for interpolation from the Poseidon-viewer (or enter them manually):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0419f74e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ps = [\n",
    "    {\"type\": \"Point\", \"coordinates\": [-57.710530333876235, 33.062070736515295]},\n",
    "    {\"type\": \"Point\", \"coordinates\": [-89.55599488497559, 25.453020491528747]},\n",
    "    {\"type\": \"Point\", \"coordinates\": [-74.70605786848537, 13.952759739624767]},\n",
    "    {\"type\": \"Point\", \"coordinates\": [-59.06355926811892, 21.999755420813273]},\n",
    "    {\"type\": \"Point\", \"coordinates\": [-41.69588507736159, 22.43955659234939]},\n",
    "    {\"type\": \"Point\", \"coordinates\": [-45.113026771172656, 14.563543828761837]},\n",
    "    {\"type\": \"Point\", \"coordinates\": [-26.081872200732885, 6.414099524482438]},\n",
    "    {\"type\": \"Point\", \"coordinates\": [-17.396102656758963, -4.381322875209875]},\n",
    "    {\"type\": \"Point\", \"coordinates\": [-26.702603318403906, -7.125636489486197]},\n",
    "    {\"type\": \"Point\", \"coordinates\": [-32.51011235240231, -22.847802807885373]},\n",
    "]\n",
    "lons, lats = ospy.utils.viewer_to_range(ps)\n",
    "Nstations = len(lons)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "046ca74e",
   "metadata": {},
   "source": [
    "## `seaduck.OceInterp` interpolates the temperature and salinity to the specified coordinates.\n",
    "\n",
    "This process makes synthetic hydrographic profiles from ECCO. Compute the potential density anomaly from the T/S data too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2df2bf6b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tub = sd.OceData(od._ds)\n",
    "times = sd.utils.convert_time(this_time)\n",
    "Ntimes = times.size\n",
    "\n",
    "depths = tub.Z.ravel()\n",
    "Ndepths = len(depths)\n",
    "\n",
    "sd_lons, sd_times, sd_depths = np.meshgrid(lons, times, depths, indexing=\"ij\")\n",
    "sd_lats, _, _ = np.meshgrid(lats, times, depths, indexing=\"ij\")\n",
    "sd_lons = np.ravel(sd_lons)\n",
    "sd_lats = np.ravel(sd_lats)\n",
    "sd_times = np.ravel(sd_times)\n",
    "sd_depths = np.ravel(sd_depths)\n",
    "\n",
    "[sd_Temp, sd_S] = np.array(\n",
    "    sd.OceInterp(tub, varList, sd_lons, sd_lats, sd_depths, sd_times)\n",
    ")\n",
    "sd_depths2 = sd_depths.reshape(Nstations, Ntimes, Ndepths).squeeze()\n",
    "sd_S = sd_S.reshape(Nstations, Ntimes, Ndepths).squeeze()\n",
    "sd_Temp = sd_Temp.reshape(Nstations, Ntimes, Ndepths).squeeze()\n",
    "\n",
    "sd_sigma0 = ospy.utils.densjmd95(sd_S, sd_Temp, 0) - 1000.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2fcebf7",
   "metadata": {},
   "source": [
    "## Plot station locations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30e4d905",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cut_od = od.subsample.cutout(\n",
    "    XRange=[min(lons) - 10, max(lons) + 10], YRange=[min(lats) - 5, max(lats) + 5]\n",
    ")\n",
    "fig = plt.figure(figsize=(12, 5))\n",
    "ax = cut_od.plot.horizontal_section(varName=\"Depth\", cmap=\"Greys_r\")\n",
    "colors = [\n",
    "    \"r\",\n",
    "    \"b\",\n",
    "    \"k\",\n",
    "    \"orange\",\n",
    "    \"darkgreen\",\n",
    "    \"m\",\n",
    "    \"yellow\",\n",
    "    \"c\",\n",
    "    \"indigo\",\n",
    "    \"magenta\",\n",
    "    \"green\",\n",
    "]\n",
    "legestations_ds = np.arange(Nstations)\n",
    "for i in range(Nstations):\n",
    "    plt.plot(\n",
    "        lons[i],\n",
    "        lats[i],\n",
    "        \"o\",\n",
    "        color=colors[i],\n",
    "        label=\"station \" + str(legestations_ds[i]),\n",
    "    )\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9d4f8b4",
   "metadata": {},
   "source": [
    "## Plot vertical T/S station profiles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92de3e3f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for i in range(Nstations):\n",
    "    plt.subplot(1, 3, 1)\n",
    "    plt.plot(\n",
    "        sd_Temp[i, :],\n",
    "        depths,\n",
    "        \"o\",\n",
    "        color=colors[i],\n",
    "        label=\"station \" + str(legestations_ds[i]),\n",
    "    )\n",
    "    plt.xlabel(\"Temperature [oC]\")\n",
    "    plt.ylabel(\"Depth [m]\")\n",
    "    plt.grid()\n",
    "\n",
    "    plt.subplot(1, 3, 2)\n",
    "    plt.plot(\n",
    "        sd_S[i, :],\n",
    "        depths,\n",
    "        \"o\",\n",
    "        color=colors[i],\n",
    "        label=\"station \" + str(legestations_ds[i]),\n",
    "    )\n",
    "    plt.xlabel(\"Salinity [g/kg]\")\n",
    "    plt.grid()\n",
    "\n",
    "    plt.subplot(1, 3, 3)\n",
    "    plt.plot(\n",
    "        sd_sigma0[i, :],\n",
    "        depths,\n",
    "        \"o\",\n",
    "        color=colors[i],\n",
    "        label=\"station \" + str(legestations_ds[i]),\n",
    "    )\n",
    "    plt.xlabel(\"Density anomaly [kg/m^3]\")\n",
    "    plt.grid()\n",
    "\n",
    "plt.subplot(1, 3, 1)\n",
    "plt.xlabel(\"Temperature [oC]\")\n",
    "plt.ylabel(\"Depth [m]\")\n",
    "plt.grid()\n",
    "\n",
    "plt.subplot(1, 3, 2)\n",
    "plt.xlabel(\"Salinity [g/kg]\")\n",
    "plt.grid()\n",
    "\n",
    "plt.subplot(1, 3, 3)\n",
    "plt.xlabel(\"Density anomaly [kg/m^3]\")\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a691241",
   "metadata": {},
   "source": [
    "## Plot T/S diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a19c06ae",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for i in range(Nstations):\n",
    "    plt.plot(\n",
    "        sd_S[i, :],\n",
    "        sd_Temp[i, :],\n",
    "        \"o\",\n",
    "        color=colors[i],\n",
    "        label=\"station \" + str(legestations_ds[i]),\n",
    "    )\n",
    "plt.xlabel(\"Salinity [g/kg]\")\n",
    "plt.ylabel(\"Temperature [oC]\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "562cecf5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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"
  },
  "toc-showmarkdowntxt": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}