book/brazilian-fires-gridded-plot-carbon-monoxide-column.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Brazilian Fires - Carbon Monoxide (CO) - CrIS JPSS-1\n",
    "\n",
    "## Total Column - gridded data (14km x 14km) plots\n",
    "\n",
    "This code plots the Carbon Monoxide (CO) total column for each day from 01-Aug-2022 to 31-Aug-2022."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "import contextlib\n",
    "from pathlib import Path\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import numpy as np\n",
    "from netCDF4 import Dataset\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "\n",
    "from scipy.interpolate import griddata\n",
    "from scipy.spatial import cKDTree\n",
    "\n",
    "import imageio.v2 as imageio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify date range\n",
    "start_date = datetime.date(2022, 8, 1)\n",
    "end_date = datetime.date(2022, 8, 31)\n",
    "\n",
    "# Loop through the date range and plot Carbon Monoxide for each day\n",
    "for current_date in pd.date_range(start_date, end_date):\n",
    "    data_file_date = current_date.strftime(\"%Y%m%d\")\n",
    "    \n",
    "    # We may not have data for certain days    \n",
    "    data_file_path = Path(f'./data/TROPESS_CrIS-JPSS1_L2_Summary_CO_{data_file_date}_MUSES_R1p17_FS_F0p6.nc')\n",
    "    if not data_file_path.exists():\n",
    "        continue\n",
    "\n",
    "    # Open the netCDF file\n",
    "    dataset = Dataset(data_file_path.as_posix(), 'r')\n",
    "\n",
    "    # Read the data from your variables\n",
    "    latitude = dataset.variables['latitude'][:]\n",
    "    longitude = dataset.variables['longitude'][:]\n",
    "    x_col = dataset.variables['x_col'][:]\n",
    "\n",
    "    dataset.close()\n",
    "\n",
    "    # Filter the data for South America only\n",
    "    latitude_max = 14.0\n",
    "    latitude_min = -26.0\n",
    "    longitude_max = -30.0\n",
    "    longitude_min = -86.0\n",
    "\n",
    "    south_america_only = np.where(\n",
    "        (latitude > latitude_min) & (latitude < latitude_max) & \n",
    "        (longitude > longitude_min) &(longitude < longitude_max)\n",
    "    )[0]\n",
    "\n",
    "    latitude = latitude[south_america_only]\n",
    "    longitude = longitude[south_america_only]\n",
    "    x_col = x_col[south_america_only]\n",
    "\n",
    "    # Calculate width and height\n",
    "    aspect_ratio = (longitude_max - longitude_min) / (latitude_max - latitude_min)\n",
    "    w = 3600; h = w / aspect_ratio\n",
    "\n",
    "    # Specify figure size (in inches)\n",
    "    dpi = 300; \n",
    "    plt.figure(figsize=(w / dpi, h / dpi))  \n",
    "\n",
    "    # Create a basemap instance\n",
    "    m = Basemap(projection='cyl', resolution='i',\n",
    "                llcrnrlat=latitude_min, urcrnrlat=latitude_max,  # set latitude limits to -56 and 24\n",
    "                llcrnrlon=longitude_min, urcrnrlon=longitude_max)  # set longitude limits to -96 and -30\n",
    "\n",
    "    m.drawmapboundary(fill_color='black')\n",
    "    m.drawcoastlines(linewidth=0.3, color='white')\n",
    "    m.drawcountries(linewidth=0.2, color='white')\n",
    "\n",
    "    # Draw parallels (latitude lines) and meridians (longitude lines)\n",
    "    parallels = np.arange(-90., 91., 10.)\n",
    "    m.drawparallels(parallels, labels=[True,False,False,False], linewidth=0.3)\n",
    "\n",
    "    meridians = np.arange(-180., 181., 10.)\n",
    "    m.drawmeridians(meridians, labels=[False,False,False,True], linewidth=0.3)\n",
    "\n",
    "    # Get CrIS pixel size in degrees, pixel is 14 x 14 km, 1 degree is roughly 111.1 km\n",
    "    cris_pixel_size_deg = 14.0 / 111.1 \n",
    "\n",
    "    # Get the grid for the interpolated values\n",
    "    grid_lat, grid_lon = np.mgrid[-26:13:cris_pixel_size_deg, -86:-29:cris_pixel_size_deg]\n",
    "\n",
    "    # Interpolate the data using griddata\n",
    "    grid_x_col = griddata((latitude, longitude), x_col, (grid_lat, grid_lon), method='linear', rescale=True)\n",
    "\n",
    "    # Find the distance to the nearest original point for each point in the interpolated grid\n",
    "    tree = cKDTree(np.vstack((latitude, longitude)).T)\n",
    "    dist, _ = tree.query(np.vstack((grid_lat.ravel(), grid_lon.ravel())).T)\n",
    "\n",
    "    # Reshape the distance array to have the same shape as the x_col grid\n",
    "    dist_grid = dist.reshape(grid_x_col.shape)\n",
    "\n",
    "    # Mask the interpolated values that are too far from any original point\n",
    "    max_distance_degrees = 3.0\n",
    "    grid_x_col[dist_grid > max_distance_degrees] = np.nan\n",
    "\n",
    "    # Plot the interpolated data using pcolormesh instead of scatter\n",
    "    sc = m.pcolormesh(grid_lon, grid_lat, grid_x_col, latlon=True, cmap='jet', alpha=0.7, vmin=60.0, vmax=200.0)\n",
    "\n",
    "    # Add a colorbar\n",
    "    cbar = m.colorbar(sc, location='bottom', pad=\"10%\")\n",
    "    cbar.set_label('Column-averaged dry-air mixing ratio of carbon monoxide from the surface to Top of Atmosphere (TOA), ppbv')\n",
    "\n",
    "    # set plot title\n",
    "    title_date = current_date.strftime(\"%d-%b-%Y\")\n",
    "    title_min = np.min(x_col)\n",
    "    title_max = np.max(x_col)\n",
    "    plt.title(f'TROPESS - CrIS JPSS-1 - Carbon Monoxide Column (XCO) - {title_date}\\nMin: {title_min:.01f}, Max: {title_max:.01f} ppbv')\n",
    "\n",
    "    # Save figure to PNG file\n",
    "    plt.savefig(f'./images/TROPESS_Brazilian-Fires_CrIS-JPSS-1_XCO_{title_date}.png', dpi=dpi, bbox_inches='tight')\n",
    "\n",
    "    plt.close()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now generate an animated GIF from the saved figures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timelapse_start_date = start_date.strftime(\"%d-%b-%Y\")\n",
    "timelapse_end_date = end_date.strftime(\"%d-%b-%Y\")\n",
    "timelapse_file = f'./images/TROPESS_Brazilian-Fires_CrIS-JPSS-1_XCO_Timelapse_{timelapse_start_date}_{timelapse_end_date}.gif'\n",
    "\n",
    "# Path to the directory where the images were saved in the previous step\n",
    "img_dir = Path('./images')\n",
    "\n",
    "# Get all file names in the directory\n",
    "files = sorted([img for img in img_dir.glob(\"TROPESS_Brazilian-Fires_CrIS-JPSS-1_XCO_*.png\")])\n",
    "\n",
    "# Read the images and add to a list\n",
    "with contextlib.redirect_stdout(None):\n",
    "    images = []\n",
    "    for filename in files:\n",
    "        images.append(imageio.imread(filename))\n",
    "\n",
    "    # Save the images as an animated gif\n",
    "    # duration is the time spent on each image (in milliseconds)\n",
    "    imageio.mimsave(timelapse_file, images, duration=1000.0, loop=0)  "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First day \n",
    "\n",
    "![first-day](./images/TROPESS_Brazilian-Fires_CrIS-JPSS-1_XCO_01-Aug-2022.png)\n",
    "\n",
    "## Timelapse\n",
    "\n",
    "![timelapse](./images/TROPESS_Brazilian-Fires_CrIS-JPSS-1_XCO_Timelapse_01-Aug-2022_31-Aug-2022.gif)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "b7cfc0ec56b1d6945f565255d60f992e290460feecef750e75b6663be732b9eb"
   }
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}