MPAS native output support + example

.gitignore

@@ -1,6 +1,7 @@
 *.ipynb_checkpoints
 Satellite_data.nc
 MPAS_data.nc
+MPAS_unstructured_data.nc
 docs/api/
 docs/examples/*.pdf
 docs/examples/*.txt

README.md

@@ -1,5 +1,5 @@
 <h1>
-  <!-- TODO: change src to https://raw.githubusercontent.com/knubez/TAMS/main/docs/_static/TAMS-logo.png once public -->
+  <!-- <img src="https://raw.githubusercontent.com/knubez/TAMS/main/docs/_static/TAMS-logo.png" -->
   <img src="./docs/_static/TAMS-logo.png"
        alt="TAMS logo" height="180" valign="bottom">
 </h1>
@@ -14,10 +14,6 @@ TAMS
 (**T**racking **A**lgorithm for [**M**esoscale Convective **S**ystems](https://en.wikipedia.org/wiki/Mesoscale_convective_system))
 in Python and with more flexibility.
 
-- Sample satellite data to test TAMS v2.0 is located at: https://drive.google.com/file/d/1HAhAlfqZGjnTk8NAjyx_lmVumUu_1TMp/view?usp=sharing
-
-- Sample MPAS data to test TAMS v2.0 is located at: https://drive.google.com/file/d/1vtx6UeSS8FM5Hy9DEQe3x78Ey-Hn-83E/view?usp=sharing
-
 - TAMS paper to cite: Núñez Ocasio et al. 2020, doi:[10.1175/MWR-D-19-0070.1](https://doi.org/10.1175/MWR-D-19-0070.1)
 
-- If you are interested in the original TAMS version (TAMS v1.0), please email knocasio@ucar.edu
+- If you are interested in the original TAMS version (TAMS v1.0), please email `knocasio@ucar.edu`

docs/api.rst

@@ -42,8 +42,9 @@ Data
 .. autosummary::
    :toctree: api/
 
-   tams.load_example_mpas
    tams.load_example_tb
+   tams.load_example_mpas
+   tams.load_example_mpas_ug
    tams.load_mpas_precip
    tams.data.download_examples
    tams.data.load_example_ir

docs/conf.py

@@ -19,6 +19,7 @@
     "sphinxext.opengraph",
     "sphinx-prompt",
     "sphinx_copybutton",
+    "sphinx_togglebutton",
 ]
 
 exclude_patterns = [

docs/examples/sample-mpas-ug-data.ipynb

@@ -0,0 +1,290 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "d7b6b1ad-3ab4-494b-a76f-93b391a3f028",
+   "metadata": {},
+   "source": [
+    "# MPAS unstructured grid data\n",
+    "\n",
+    "In this example, we demonstrate CE identification with {func}`unstructured-grid data <tams.load_example_mpas_ug>` from MPAS\n",
+    "and compare to the results using the {func}`regridded dataset <tams.load_example_mpas>`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9c3986a8-51f9-4a8e-bdd7-01f6672a780c",
+   "metadata": {
+    "tags": [
+     "hide-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "import warnings\n",
+    "\n",
+    "import cartopy\n",
+    "import cartopy.crs as ccrs\n",
+    "import cartopy.feature as cfeature\n",
+    "import matplotlib as mpl\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import xarray as xr\n",
+    "from matplotlib.tri import Triangulation\n",
+    "\n",
+    "import tams\n",
+    "\n",
+    "warnings.filterwarnings(\"ignore\", category=cartopy.io.DownloadWarning)\n",
+    "\n",
+    "%matplotlib inline\n",
+    "\n",
+    "xr.set_options(display_expand_data=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3bb91273-9857-479a-b2f4-be502cc52c1a",
+   "metadata": {},
+   "source": [
+    "## Load data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d552c7c5-b5fe-4bdb-b324-738a29453b78",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = tams.load_example_mpas_ug()\n",
+    "ds"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5483b3a7-025e-47e7-9ff2-1c306fce00e8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(6, 3))\n",
+    "\n",
+    "sel = ds.isel(cell=slice(None, None, 20))\n",
+    "ax.scatter(sel.lon, sel.lat, marker=\".\", s=10, alpha=0.5, edgecolors=\"none\")\n",
+    "ax.set(xlabel=\"lon\", ylabel=\"lat\")\n",
+    "ax.autoscale(tight=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e416da1f-ee50-4887-afcf-c79a552f279c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(7.5, 3))\n",
+    "\n",
+    "im = ax.scatter(ds.lon, ds.lat, c=ds.tb.isel(time=10), marker=\".\", s=3, edgecolors=\"none\")\n",
+    "fig.colorbar(im, ax=ax, label=\"Tb\")\n",
+    "ax.autoscale(tight=True)\n",
+    "ax.set(xlabel=\"lon\", ylabel=\"lat\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6b852b14-b64d-4d20-9177-a224938ffab4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ts = ds.mean(\"cell\", keep_attrs=True)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(6, 3))\n",
+    "ax2 = ax.twinx()\n",
+    "\n",
+    "ts.tb.plot(ax=ax, c=\"orangered\")\n",
+    "ts.precip.plot(ax=ax2, c=\"cornflowerblue\")\n",
+    "ax.autoscale(axis=\"x\", tight=True)\n",
+    "ax.grid(True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4b310c93-5c28-49aa-ab9c-2da1a3b66fe8",
+   "metadata": {},
+   "source": [
+    "## Identify CEs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e6df13ec-8401-4613-9368-a0aad169522c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "\n",
+    "itime = 10\n",
+    "\n",
+    "stime = pd.Timestamp(ds.time.values[itime]).strftime(r\"%Y-%m-%d_%H\")\n",
+    "print(stime)\n",
+    "\n",
+    "x = ds.lon\n",
+    "y = ds.lat\n",
+    "tri = Triangulation(x, y)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3b78cb23-6cd4-4838-87c4-2a7b31780154",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "\n",
+    "# Passing the triangulation in is not required but makes it faster\n",
+    "shapes = tams.core._contours_to_gdf(tams.contours(ds.tb.isel(time=itime), value=235, triangulation=tri));"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "87cd75b7-7563-4b74-8ebb-d4cf19c6e953",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "\n",
+    "# `tams.identify` does the above but also for the core threshold and does size filtering by default\n",
+    "cs_ug, cs_ug_core = list(zip(*tams.identify(ds.tb.isel(time=10))))[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "95225bb2-0706-421c-bb92-1478f095b509",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "\n",
+    "tran = ccrs.PlateCarree()\n",
+    "proj = ccrs.PlateCarree()  # near equator\n",
+    "\n",
+    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6), subplot_kw=dict(projection=proj), constrained_layout=True)\n",
+    "\n",
+    "ax = ax1\n",
+    "\n",
+    "ax.add_feature(cfeature.LAND)\n",
+    "\n",
+    "tcs = ax.tricontour(tri, ds.tb.isel(time=itime), levels=[235], colors=\"red\", linewidths=1, transform=tran)\n",
+    "tcs_core = ax.tricontour(tri, ds.tb.isel(time=itime), levels=[219], colors=\"blue\", linewidths=1, transform=tran)\n",
+    "\n",
+    "shapes.plot(fc=\"none\", ec=\"green\", lw=1.5, ls=\":\", ax=ax, transform=tran)  # not size-filtered\n",
+    "cs_ug.plot(fc=\"none\", ec=\"lawngreen\", lw=2, ax=ax, zorder=3, transform=tran)\n",
+    "\n",
+    "ax.gridlines(draw_labels=True)\n",
+    "\n",
+    "legend_handles = [\n",
+    "    mpl.patches.Patch(color=\"red\", label=\"219 K contours\"),\n",
+    "    mpl.patches.Patch(color=\"blue\", label=\"235 K contours\"),\n",
+    "    mpl.patches.Patch(color=\"lawngreen\", label=\"CE polygons\"),\n",
+    "]\n",
+    "ax.legend(handles=legend_handles, loc=\"upper right\")\n",
+    "\n",
+    "ax = ax2\n",
+    "\n",
+    "ds_rg = tams.load_example_mpas().sel(lat=slice(None, 20))  # same lat upper bound as in the ug data\n",
+    "\n",
+    "ax.add_feature(cfeature.LAND)\n",
+    "\n",
+    "cs_rg, cs_rg_core = list(zip(*tams.identify(ds_rg.tb.isel(time=itime))))[0]\n",
+    "\n",
+    "a = cs_rg_core.plot(fc=\"none\", ec=\"royalblue\", lw=1.5, ls=\"--\", transform=tran, ax=ax)\n",
+    "cs_ug_core.plot(fc=\"none\", ec=\"mediumblue\", lw=2, transform=tran, ax=ax)\n",
+    "\n",
+    "cs_rg.plot(fc=\"none\", ec=\"red\", lw=1.5, ls=\"--\", transform=tran, ax=ax)\n",
+    "cs_ug.plot(fc=\"none\", ec=\"firebrick\", lw=2, zorder=3, transform=tran, ax=ax)\n",
+    "\n",
+    "legend_handles = [\n",
+    "    mpl.lines.Line2D([], [], color=\"royalblue\", ls=\"--\", lw=1.5, label=\"Regridded | cold-core polygons\"),\n",
+    "    mpl.lines.Line2D([], [], color=\"red\", ls=\"--\", lw=1.5, label=\"Regridded | CE polygons\"),\n",
+    "    mpl.lines.Line2D([], [], color=\"mediumblue\", ls=\"-\", lw=2, label=\"Native | cold-core polygons\"),\n",
+    "    mpl.lines.Line2D([], [], color=\"firebrick\", ls=\"-\", lw=2, label=\"Native | CE polygons\"),\n",
+    "]\n",
+    "\n",
+    "gl = ax.gridlines(draw_labels=True)\n",
+    "gl.top_labels = False\n",
+    "# ax.autoscale(tight=True)\n",
+    "ax.set_xlim(ax1.get_xlim()); ax.set_ylim(ax1.get_ylim())\n",
+    "fig.legend(handles=legend_handles, ncol=2, loc=\"lower right\", bbox_to_anchor=[0, -0.09, 0.961, 1], frameon=False)\n",
+    "\n",
+    "for a, ax in zip(\"abc\", [ax1, ax2]):\n",
+    "    ax.text(0.007, 0.98, a, weight=\"bold\", size=14, va=\"top\", ha=\"left\", transform=ax.transAxes)\n",
+    "\n",
+    "fig.savefig(f\"mpas-ug-contours-and-vs-rg-ces_{stime}.pdf\", bbox_inches=\"tight\", pad_inches=0.05, transparent=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "df1e4cfe",
+   "metadata": {},
+   "source": [
+    "## Triangulations\n",
+    "\n",
+    "The above contouring is based on the [Delauney triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation).\n",
+    "This connects the cell centers with straight lines to form triangles.\n",
+    "But this triangulation doesn't correspond to the MPAS mesh, per se.\n",
+    "\n",
+    "MPAS-Tools, though, [provides a method](https://mpas-dev.github.io/MPAS-Tools/stable/visualization.html#mpas-mesh-to-triangles) for dividing the mesh grid cells into triangles\n",
+    "and interpolating data defined at cell centers (like OLR) to the triangle nodes.\n",
+    "The triangle nodes include the cell centers (one per triangle)\n",
+    "and the grid cell vertices (two per triangle).\n",
+    "A hexagonal grid cell, e.g., is divided into 6 triangles,\n",
+    "which all share the cell center as a node.\n",
+    "\n",
+    "We compared the contourings from these two methods.\n",
+    "For our intents and purposes\n",
+    "(identifying CEs with area ≥ 4000 km²),\n",
+    "the differences appear to be negligible.\n",
+    "The mesh horizontal resolution is 15 km,\n",
+    "which seems to be sufficiently high that the Delauney triangulation is a good approximation.\n",
+    "\n",
+    "```{code-block} python\n",
+    ":class: toggle\n",
+    "\n",
+    "# Code\n",
+    "```\n",
+    "\n",
+    "![tricontour comparison zoomed](../_static/tri-compare_zoom.png)"
+   ]
+  }
+ ],
+ "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",
+   "version": "3.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/examples/sample-satellite-data.ipynb

@@ -6,13 +6,17 @@
    "source": [
     "# Sample satellite data\n",
     "\n",
-    "Some graphical tests of current primary functions using the sample satellite data."
+    "Some graphical tests of current primary functions using the {func}`sample satellite data <tams.load_example_tb>`."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "tags": [
+     "hide-output"
+    ]
+   },
    "outputs": [],
    "source": [
     "import warnings\n",
@@ -24,16 +28,9 @@
     "\n",
     "import tams\n",
     "\n",
-    "warnings.filterwarnings(\"ignore\", category=cartopy.io.DownloadWarning)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.close(\"all\")\n",
+    "warnings.filterwarnings(\"ignore\", category=cartopy.io.DownloadWarning)\n",
+    "\n",
+    "%matplotlib inline\n",
     "\n",
     "xr.set_options(display_expand_data=False)"
    ]
@@ -61,7 +58,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "tb.isel(time=0).plot(x=\"lon\", y=\"lat\", size=5, aspect=2.5)"
+    "tb.isel(time=0).plot(x=\"lon\", y=\"lat\", size=5, aspect=2.5);"
    ]
   },
   {
@@ -102,7 +99,7 @@
     "fig, ax = plt.subplots()\n",
     "c = contour_sets[m].iloc[[n]]\n",
     "c.plot(ax=ax)\n",
-    "c.cs219.plot(color=\"red\", ax=ax, alpha=0.4)"
+    "c.cs219.plot(color=\"red\", ax=ax, alpha=0.4);"
    ]
   },
   {