hurricane-florence-animation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d9ff4103",
   "metadata": {},
   "source": [
    "## Visualizing Hurricane Florence\n",
    "\n",
    "This examples makes a true-color animation of Hurricane Florence by stitching together images from GOES. It builds off this example from [pytroll-examples](https://github.com/pytroll/pytroll-examples/blob/main/satpy/GOES-16%20ABI%20-%20True%20Color%20Animation%20-%20Hurricane%20Florence.ipynb). You can see the output of that example [here](https://twitter.com/PyTrollOrg/status/1039555399433834497).\n",
    "\n",
    "Here's what our final animation will look like:\n",
    "\n",
    "<video controls width=\"80%\" src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/pc_video/pc-examples-goes-florence.webm\">\n",
    "</video>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d48a3b9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import urllib.request\n",
    "\n",
    "import contextily as ctx\n",
    "import geopandas\n",
    "import matplotlib.animation as animation\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import planetary_computer\n",
    "import pystac_client\n",
    "import rioxarray\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6217412c",
   "metadata": {},
   "source": [
    "### Find the storm\n",
    "\n",
    "First, we need to find where on earth the storm was. The NCEI [International Best Track Archive for Climate Stewardship](https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.ncdc:C00834) provides files with all the information we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9601ba0",
   "metadata": {},
   "outputs": [],
   "source": [
    "file, _ = urllib.request.urlretrieve(\n",
    "    \"https://www.ncei.noaa.gov/data/international-best-track-archive-for-\"\n",
    "    \"climate-stewardship-ibtracs/v04r00/access/netcdf/IBTrACS.NA.v04r00.nc\"\n",
    ")\n",
    "# The storm id comes from the text file in\n",
    "# https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs\n",
    "# /v04r00/access/netcdf/\n",
    "# The name of this file changes with the update date, so we can't access it programatically.\n",
    "STORM_ID = b\"2018242N13343\"\n",
    "ds = xr.open_dataset(file)\n",
    "storm_loc = (ds.sid == STORM_ID).argmax().item()\n",
    "\n",
    "data = ds.sel(storm=storm_loc)\n",
    "geometry = geopandas.points_from_xy(data.lon, data.lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5b46235",
   "metadata": {},
   "source": [
    "`geometry` is a geopandas GeoArray with points tracking the location of the storm over time. We'll match those up with the timestamps to plot plot storm's trajectory. We'll also overlay the time period covered by our animation in red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b5ada5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = (\n",
    "    geopandas.GeoDataFrame(\n",
    "        dict(\n",
    "            time=pd.to_datetime(data.time).tz_localize(\"UTC\"),\n",
    "            geometry=geopandas.points_from_xy(data.lon, data.lat),\n",
    "        )\n",
    "    )\n",
    "    .set_crs(4326)\n",
    "    .dropna()\n",
    ")\n",
    "\n",
    "start = pd.Timestamp(\"2018-09-11T13:00:00Z\")\n",
    "stop = pd.Timestamp(\"2018-09-11T15:40:00Z\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbd23142",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = df.to_crs(epsg=3857).plot(figsize=(12, 12))\n",
    "subset = df[df.time.dt.date == start.date()]\n",
    "subset.to_crs(epsg=3857).plot(ax=ax, color=\"r\")\n",
    "\n",
    "ctx.add_basemap(ax)\n",
    "ax.set_axis_off()\n",
    "ax.set(title=\"Path of Hurricane Florence (animation period in red)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ca49d78",
   "metadata": {},
   "source": [
    "Let's save the bounding box for the subset of points we're animating . We'll use it in our query later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc950ea3",
   "metadata": {},
   "outputs": [],
   "source": [
    "bbox = list(subset.total_bounds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1330e644",
   "metadata": {},
   "source": [
    "### Get the imagery\n",
    "\n",
    "Now we'll get the GOES imagery using the Planteary Computer's STAC API. We'll use the `goes-cmi` collection. We'll also have the API filter down the images to just the \"mesoscale\" images (GOES takes images with various fields of view)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4df7dad0",
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog = pystac_client.Client.open(\n",
    "    \"https://planetarycomputer.microsoft.com/api/stac/v1/\"\n",
    ")\n",
    "search = catalog.search(\n",
    "    collections=[\"goes-cmi\"],\n",
    "    bbox=bbox,\n",
    "    datetime=[start, stop],\n",
    "    limit=500,\n",
    "    query={\"goes:image-type\": {\"eq\": \"MESOSCALE\"}},\n",
    ")\n",
    "items = search.get_all_items()\n",
    "signed_items = sorted(\n",
    "    [planetary_computer.sign(item) for item in items], key=lambda x: x.datetime\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "333e4d59",
   "metadata": {},
   "source": [
    "Let's load and plot the first item, just to make sure we're on the right track."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92855bfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = rioxarray.open_rasterio(signed_items[0].assets[\"C01_2km\"].href).load()\n",
    "ds[0].plot.imshow(figsize=(16, 9), cmap=\"Blues\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "377a0ab0",
   "metadata": {},
   "source": [
    "Great. Now we can load all the data for the first three bands (blue, red, and near-infrared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82b0addf",
   "metadata": {},
   "outputs": [],
   "source": [
    "bands = [\"C01_2km\", \"C02_2km\", \"C03_2km\"]\n",
    "common_names = [\n",
    "    items[0].assets[band].extra_fields[\"eo:bands\"][0][\"common_name\"] for band in bands\n",
    "]\n",
    "time = xr.DataArray(\n",
    "    pd.to_datetime([x.datetime for x in signed_items]).tz_localize(None),\n",
    "    name=\"time\",\n",
    "    dims=[\"time\"],\n",
    ")\n",
    "arrays = [\n",
    "    xr.concat(\n",
    "        [rioxarray.open_rasterio(item.assets[band].href) for band in bands], dim=\"band\"\n",
    "    ).assign_coords(band=common_names)\n",
    "    for item in signed_items\n",
    "]\n",
    "data = xr.concat(arrays, dim=time).rename(\"goes\")\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3f1a578",
   "metadata": {},
   "source": [
    "Goes doesn't have a true green band, which we need for our true color animation. We'll simulate it with a linear combination of the other bands (See [Bah et. al (2018)](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2018EA000379) for more on this technique)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "733c325f",
   "metadata": {},
   "outputs": [],
   "source": [
    "green = (\n",
    "    0.45 * data.sel(band=\"red\")\n",
    "    + 0.1 * data.sel(band=\"nir09\")\n",
    "    + 0.45 * data.sel(band=\"blue\")\n",
    ").assign_coords(band=\"green\")\n",
    "green"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "107c4dba",
   "metadata": {},
   "source": [
    "Now we'll normalize the data and apply a gamma correction for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2a0a541",
   "metadata": {},
   "outputs": [],
   "source": [
    "γ = 2.2\n",
    "\n",
    "rgb = xr.concat([data, green], dim=\"band\").sel(band=[\"red\", \"green\", \"blue\"])\n",
    "rgb = rgb / rgb.max(dim=[\"band\", \"y\", \"x\"])\n",
    "rgb = np.clip(rgb ** (1 / γ), 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbb5b41b",
   "metadata": {},
   "source": [
    "Let's check out the first image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4b9fee9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 16))\n",
    "rgb.isel(time=0).plot.imshow(rgb=\"band\", add_labels=False)\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "354e41ba",
   "metadata": {},
   "source": [
    "### Create the animation\n",
    "\n",
    "We'll use matplotlib's [FuncAnimation](https://matplotlib.org/stable/api/_as_gen/matplotlib.animation.FuncAnimation.html) to create the animation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1e2adaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 16))\n",
    "ax.set_axis_off()\n",
    "\n",
    "img = rgb[0].plot.imshow(ax=ax, add_colorbar=False, rgb=\"band\", add_labels=False)\n",
    "label = ax.text(\n",
    "    0.4,\n",
    "    0.03,\n",
    "    pd.Timestamp(rgb.time.data[0]).isoformat(),\n",
    "    transform=ax.transAxes,\n",
    "    color=\"k\",\n",
    "    size=20,\n",
    ")\n",
    "\n",
    "\n",
    "def animate(i):\n",
    "    img.set_data(rgb[i].transpose(\"y\", \"x\", \"band\"))\n",
    "    label.set_text(pd.Timestamp(rgb.time.data[i]).isoformat())\n",
    "    return img, label\n",
    "\n",
    "\n",
    "ani = animation.FuncAnimation(fig, animate, frames=len(rgb), interval=120)\n",
    "ani.save(\"goes.mp4\", fps=15, extra_args=[\"-vcodec\", \"libx264\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f38385c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Video\n",
    "\n",
    "Video(\"goes.mp4\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "593fad1a",
   "metadata": {},
   "source": [
    "### Next steps\n",
    "\n",
    "Learn more about GOES and using the Planetary Computer\n",
    "\n",
    "* [GOES quickstart](../datasets/goes/goes-example.ipynb)\n",
    "* [Reading from the STAC API](https://planetarycomputer.microsoft.com/docs/quickstarts/reading-stac/)"
   ]
  }
 ],
 "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}