{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "lm9azastWMN7" }, "source": [ "# Exploring TC PRIMED, Chapter 4: Regridding Satellite Swath Data\n", "- Creators: Kathy Haynes, Naufal Razin, and Chris Slocum\n", "- Affiliations: CIRA and NESDIS/STAR\n", "\n", "___\n", "\n", "## Overview\n", "The Tropical Cyclone Precipitation, Infrared, Microwave, and Environmental Dataset (TC PRIMED) contains satellite passive microwave data in the original swath format. To apply machine learning approaches to swath data, we must either apply it on a pixel-by-pixel basis or pre-process the swath data. If you've been following along with the notebooks in this series, you would have learned to read, visualize, and analyze TC PRIMED data, as well as how to utilize it for a pixel-based machine learning application. In this notebook, we explore one way to pre-process TC PRIMED swath data for image-based machine learning applications by regridding swath data onto a uniform grid.\n", "\n", "## Prerequisites\n", "To successfully navigate and use this notebook, you should be accustomed with:\n", "- Basics of Python programming (e.g., loading modules, variable assignment, plotting, familiarity with scientific Python packages)\n", "- Basics of satellite data (e.g., familiarity with satellite swath data such as that from [low-Earth-orbiting](https://www.esa.int/Enabling_Support/Space_Transportation/Types_of_orbits) satellites)\n", "- [Map projections](https://gisgeography.com/map-projections/) and grid types (e.g., Mercator, Lambert Equal Area)\n", "- [Interpolation methods](https://gisgeography.com/raster-resampling/) (e.g., Nearest Neighbor, Gaussian Weighting, Bilinear)\n", "\n", "## Learning Outcomes\n", "After working through this notebook, you should be able to:\n", "- understanding basic satellite swath data (such as that from polar-orbiting satellites)\n", "- execute a remapping of satellite swath data onto a discrete grid" ] }, { "cell_type": "markdown", "metadata": { "id": "61MOjwekWMN-" }, "source": [ "## Background\n", "\n", "There are many types of satellite data in meteorological use.\n", " Geostationary satellites have onboard sensors that view the same portion of the Earth's surface at all times. In contrast, microwave sensors are found on polar orbiting satellites. With polar orbiters, the satellite passes over certain portions of the Earth's surface as it travels around the Earth, poleward on one side and then equatorward on the opposite side, creating 'stripes' or swaths of data. During the satellite's travels, the Earth rotates underneath it, allowing the satellite swaths to cover new areas with each consecutive pass.\n", "\n", "TC PRIMED is a dataset centered around satellite passive microwave observations of tropical cyclones, and these data are in the form of swaths. The nature of swath data adds complexity for machine learning applications using TC PRIMED, particularly for predicting storm-scale phenomena like the presence of an eye or the current intensity. Because satellite swath data only covers a portion of the Earth at one time, putting it on a regular grid allows the data to be in a standard, consistent format for machine learning. Additionally, it allows for the combination and collocation of swath data with other satellite data, allowing you to use different types of satellite data together.\n", "\n", "In this notebook, we will use an existing Python package, `pyresample`, to regrid TC PRIMED swath data onto a uniform grid. This regridding process can serve as a pre-processing step in applying more complex machine learning approaches like convolutional neural networks to TC PRIMED data." ] }, { "cell_type": "markdown", "metadata": { "id": "Ol_Js_5IWMN-" }, "source": [ "### Hurricane Florence (2018)\n", "In this example, we will use TC PRIMED files from Hurricane Florence (2018). If you have gone through the earlier notebooks in this series, you should be familiar with Hurricane Florence. Hurricane Florence was a long-lived hurricane that tracked across the Atlantic to make landfall along the southeastern coast of North Carolina as a high-end Category 1 hurricane." ] }, { "cell_type": "markdown", "metadata": { "id": "-d8AsM5eWMN_" }, "source": [ "
\n", "\"An\n", "\n", "
Figure 1. An image of Hurricane Florence (2018) captured by European Space Agency astronaut Alexander Gerst from aboard the International Space Station. Source: ESA
\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "NQ1YNKNlWMOP" }, "source": [ "### Passive microwave observations in the 85 to 92 GHz frequency window\n", "In this notebook, we will regrid passive microwave observations in the 85 to 92 GHz frequency window. Observations in this microwave frequency window can penetrate through the cloud tops to reveal the internal storm structure. In essence, it \"sees\" through cloud tops. In this frequency window, land surfaces appear warm relative to water. In addition, low-level moist air masses act to warm the brightness temperatures over water surfaces, which highlights the warm emission regions of the hurricane where there are rain and cloud water bands. In contrast, deep convection appears relatively cold due to scattering from ice, highlighting the cold scattering regime of the hurricane and areas of deep convection." ] }, { "cell_type": "markdown", "metadata": { "id": "vRUjwrGOWMOP" }, "source": [ "## Software\n", "This tutorial uses the Python programming language and scientific Python packages. We will primarily use:\n", "* `netCDF4` and `numpy` to read in satellite swath data\n", "* `pyresample` to grid the data\n", "* `cartopy` and `matplotlib` for plotting" ] }, { "cell_type": "markdown", "metadata": { "id": "22QQHcUBWMOP" }, "source": [ "### Install Packages\n", "\n", "Let's first check if we have the necessary Python packages to run this notebook. If we don't, let's install them." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 1000 }, "executionInfo": { "elapsed": 56986, "status": "ok", "timestamp": 1686768924391, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "0xhHH9X9td4v", "outputId": "8a42f449-c929-462d-e304-26b5c8b54854" }, "outputs": [], "source": [ "import subprocess, sys\n", "packages = [\"netCDF4\", \"numpy\", \"matplotlib\",\n", " \"cartopy\", \"pyresample\"]\n", "# If this notebook is running on Google Colab, we know what packages\n", "# need to be installed. In addition, we would have to resolve some\n", "# incompatibility issues between shapely and cartopy.\n", "if 'google.colab' in sys.modules:\n", " !pip install netCDF4\n", " !pip install --upgrade numpy\n", " !pip install --upgrade matplotlib\n", " !pip install --no-binary shapely shapely --force\n", " !pip install cartopy\n", " !pip install pyresample\n", "\n", "# If this note book is not running on Google Colab, check if Python\n", "# has the packages we need to run this notebook. If not, download it.\n", "else:\n", " for package in packages:\n", " try:\n", " __import__(package)\n", " except ImportError:\n", " subprocess.check_call([sys.executable, '-m', 'pip', 'install', package])" ] }, { "cell_type": "markdown", "metadata": { "id": "snu-Q4KTz96O" }, "source": [ "Now, let's load the modules in the packages (e.g., `DownloadWarning`) or load the packages and assign a shorter object name for the packages (e.g., `as ccrs`) for a cleaner use throughout the notebook." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "i_BQfL7GWMOQ" }, "outputs": [], "source": [ "# Install and import required libraries\n", "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeature\n", "from cartopy.io import DownloadWarning\n", "\n", "import matplotlib.pyplot as plt\n", "from matplotlib import cm, ticker\n", "from matplotlib.colors import Normalize\n", "\n", "from netCDF4 import Dataset\n", "import numpy as np\n", "\n", "from pyresample import kd_tree\n", "from pyresample import SwathDefinition\n", "from pyresample.area_config import create_area_def\n", "from pyresample.bilinear import NumpyBilinearResampler\n", "from pyresample.geometry import AreaDefinition\n", "from pyresample.plot import show_quicklook\n", "from pyresample.utils import check_and_wrap" ] }, { "cell_type": "markdown", "metadata": { "id": "qUDW6hoUhER7" }, "source": [ "Let's also just suppress any warnings so that we don't have to look through those." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "py6LEyY3gugc" }, "outputs": [], "source": [ "# Import warnings and suppress relevant categories.\n", "import warnings\n", "warnings.simplefilter(action='ignore', category=FutureWarning)\n", "warnings.simplefilter(action='ignore', category=DeprecationWarning)\n", "warnings.simplefilter(action='ignore', category=DownloadWarning)\n", "warnings.simplefilter(action='ignore', category=RuntimeWarning)\n", "warnings.simplefilter(action='ignore', category=UserWarning)" ] }, { "cell_type": "markdown", "metadata": { "id": "RFVWE06YM98K" }, "source": [ "## Read File Online\n", "Finally, let's retrieve information from the TC PRIMED file that we will use in this example. The TC PRIMED file will be from an Advanced Scanning Microwave Radiometer 2 (AMSR2) overpass of Hurricane Florence (2018). We have selected an AMSR2 overpass from Hurricane Florence with good coverage of the storm while it made landfall. We will use the Python `netCDF4` and `requests` packages to read and retrieve the information directly from the TC PRIMED file available on an Amazon Web Service S3 bucket as part of the NOAA Open Data Dissemination program (NODD), without downloading the file, and store the information from the file in an \"instance\" type called `DS`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "j3XDeeRSnpps" }, "outputs": [], "source": [ "import requests\n", "\n", "# Specify the URL to the TC PRIMED folder on NODD\n", "NODD_URL = \"https://noaa-nesdis-tcprimed-pds.s3.amazonaws.com/v01r00/final/2018/AL/06/\"\n", "\n", "# Specify the name of the file we will use from the TC PRIMED folder on NODD\n", "FILE_NAME = \"TCPRIMED_v01r00-final_AL062018_AMSR2_GCOMW1_033650_20180914065914.nc\"\n", "\n", "# Join NODD_URL and FILE_NAME to produce a complete link\n", "# Retrieve the contents of the TC PRIMED file from the complete link\n", "url_response = requests.get(NODD_URL + FILE_NAME)\n", "\n", "# Load the contents of the TC PRIMED file in an \"instance\" called DS\n", "DS = Dataset(FILE_NAME, memory=url_response.content)" ] }, { "cell_type": "markdown", "metadata": { "id": "szwMkO0feTK0" }, "source": [ "Now that the file is available, let's pull out the necessary information we'll need. Chapter 1b in this notebook series outlined in detail the structure of the TC PRIMED overpass file, so here we will quickly obtain the storm information along with the swath latitude, longitude, and the 89 GHz brightness temperatures.\n", "\n", "For the longitude, TC PRIMED uses longitudes from 0 to 360; however, Pyresample needs them to be from -180 to 180. Here we adjust these to make the longitude in the correct range for regridding." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "executionInfo": { "elapsed": 117, "status": "ok", "timestamp": 1686768951196, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "6XTG67yFfIOS", "outputId": "d4f49fc5-76ee-4d16-96c5-f7015e71f23f" }, "outputs": [], "source": [ "# Pull out the storm's center location\n", "storm_latitude = DS[\"overpass_storm_metadata/storm_latitude\"][:][0]\n", "storm_longitude = DS[\"overpass_storm_metadata/storm_longitude\"][:][0]\n", "if storm_longitude > 180:\n", " storm_longitude -= 360.\n", "print(\"Storm Center Longitude= {:.2f} E, Latitude= {:.2f} N\".\\\n", " format(storm_longitude, storm_latitude))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "zttkZkRxgoMr" }, "outputs": [], "source": [ "# Pull the swath latitude and longitude from the S6 swath\n", "swath_latitude = DS[\"passive_microwave/S6/latitude\"][:]\n", "swath_longitude = DS[\"passive_microwave/S6/longitude\"][:]\n", "swath_longitude[swath_longitude >= 180] -= 360" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "executionInfo": { "elapsed": 154, "status": "ok", "timestamp": 1686771869019, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "K77xhQUghIQm", "outputId": "3959f401-df41-4d61-a271-7d4af52fa894" }, "outputs": [], "source": [ "# Pull the horizontally-polarized 89 GHz brightness temperature\n", "# from the S6 swath\n", "swath_h89 = DS[\"passive_microwave/S6/TB_B89.0H\"][:]\n", "print(\"89-GHz Brightness Temperature Minimum= {:.2f}, Maximum= {:.2f}\".\\\n", " format(swath_h89.min(), swath_h89.max()))" ] }, { "cell_type": "markdown", "metadata": { "id": "BcN4NvqEcD-O" }, "source": [ "## Close the File\n", "When loading data from a NetCDF file, **always remember to close the file**. Since we will no longer be retrieving data from the file, we can close the file now." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "dJLB4A0MtoT_" }, "outputs": [], "source": [ "# Close the file\n", "DS.close()" ] }, { "cell_type": "markdown", "metadata": { "id": "Ke3jWaCwiHpe" }, "source": [ "## Visualize Swath Data\n", "\n", "In order to visualize the passive microwave data throughout this notebook, let's define a function to map the data onto various user-specified projections using `cartopy`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Q19-S97Ahoim" }, "outputs": [], "source": [ "# Define function to map 2D data\n", "def map2D(data2D, lon2D, lat2D,\n", " contours=None,\n", " cbar=True,\n", " cmap='viridis',\n", " cFontSize=16,\n", " cLabel='K',\n", " cPad=0.0,\n", " cTicks=None,\n", " dLat=0.,\n", " dLon=0.,\n", " extent=None,\n", " figSize=(10, 6),\n", " fontSize=24,\n", " gridLines=True,\n", " gridLineLabels=True,\n", " labelFormat=None,\n", " nContours=42,\n", " projection='Mercator',\n", " saveFile='',\n", " title='',\n", " vmax=None,\n", " vmin=None\n", " ):\n", " \"\"\"\n", " Create a map of 2D data with specified latitude and longitude locations.\n", " Produces an inline figure or a saved file (if saveFile is a file name).\n", "\n", " Parameters\n", " ----------\n", " data2D : array-like\n", " array of 2D data to visualize\n", " lon2D : array-like\n", " array of 2D longitudes for the data\n", " lat2D : array-like\n", " array of 2D latitudes for the data\n", " contours : array-like, optional\n", " array of contour levels for plot, calculated from data\n", " and nContours if None\n", " cbar : bool, default=True\n", " flag to show colorbar\n", " cmap : str, default='viridis'\n", " color table of data\n", " cFontSize : int, default=16\n", " font size of colorbar\n", " cLabel : str, default='K'\n", " colorbar label\n", " cPad : float, default=0.0\n", " padding for colorbar\n", " cTicks : array-like, optional\n", " specific tick locations for colorbar\n", " dLat : float, default=0.\n", " extra latitudes to show outside of data area (degrees)\n", " dLon : float, default=0.\n", " extra longitudes to show outside of data area (degrees)\n", " extent : array-like, optional\n", " extent of area to show, calculated from data and dLat/dLon if None\n", " figSize : tuple, default=(10, 6)\n", " figure size\n", " fontSize : int, default=24\n", " font size of figure title\n", " gridLines : bool, default=True\n", " flag to show grid lines\n", " gridLineLabels : bool, default=True\n", " flag to show grid line labels\n", " labelFormat : f-str, default=None\n", " format for colorbar tick labels\n", " nContours : int, default=42\n", " number of colorbar contours\n", " projection : str, default='Mercator'\n", " projection of background map\n", " saveFile : str, default=''\n", " name of file to save image, image saved only if not ''\n", " title : str, default=''\n", " figure title, printed only if not ''\n", " vmax : float, optional\n", " minimum value for colorbar, calculated from data if None\n", " vmin : float, optional\n", " maximum value for colorbar, calculated from data if None\n", "\n", " Yields\n", " ------\n", " image of regridded satellite swath data on a map with a projection\n", " \"\"\"\n", "\n", "\n", " # Get valid data ranges\n", " myMin = np.nanmin(data2D)\n", " myMax = np.nanmax(data2D)\n", " if contours is None:\n", " cDelta = (myMax - myMin) / (nContours - 2)\n", " contours = np.empty(nContours)\n", " for i in range(nContours):\n", " contours[i] = myMin + i * cDelta - 0.5 * cDelta\n", "\n", " # Setup map projection using cartopy\n", " cLat = lat2D.min() + (lat2D.max() - lat2D.min()) * 0.5\n", " cLon = lon2D.min() + (lon2D.max() - lon2D.min()) * 0.5\n", "\n", " if projection == 'LambertConformal':\n", " crs = ccrs.LambertConformal(\n", " central_longitude=cLon,\n", " central_latitude=cLat)\n", " elif projection == 'PlateCarree':\n", " crs = ccrs.PlateCarree(central_longitude=cLon)\n", " elif projection == 'Mollweide':\n", " crs = ccrs.Mollweide(central_longitude=cLon)\n", " elif projection == 'Robinson':\n", " crs = ccrs.Robinson(central_longitude=cLon)\n", " elif projection == 'Sinusoidal':\n", " crs = ccrs.Sinusoidal(central_longitude=cLon)\n", " elif projection == 'LambertConformal':\n", " crs = ccrs.LambertConformal()\n", " elif projection == 'AlbersEqualArea':\n", " crs = ccrs.AlbersEqualArea(\n", " central_longitude=cLon,\n", " central_latitude=cLat)\n", " elif projection == 'LambertAzimuthalEqualArea':\n", " crs = ccrs.LambertAzimuthalEqualArea(\n", " central_longitude=cLon,\n", " central_latitude=cLat)\n", " elif projection == 'Mercator':\n", " crs = ccrs.Mercator(central_longitude=cLon)\n", " else:\n", " print(\"Unknown projection, using PlateCarree.\")\n", " crs = ccrs.PlateCarree(central_longitude=cLon)\n", "\n", " # Create the figure and plot the background\n", " fig, ax = plt.subplots(ncols=1, nrows=1,\n", " figsize=figSize,\n", " subplot_kw={'projection': crs})\n", "\n", " if extent is None:\n", " extent = [lon2D.min() - dLon, lon2D.max() + dLon,\n", " lat2D.min() - dLat, lat2D.max() + dLat]\n", " ax.set_extent(extent)\n", " if gridLines:\n", " gl = ax.gridlines(draw_labels=gridLineLabels)\n", " gl.bottom_labels = True\n", " gl.left_labels = True\n", " gl.top_labels = False\n", " gl.right_labels = False\n", "\n", " # Draw the maps\n", " if labelFormat is None:\n", " labelFormat = '%.0f'\n", " cf = ax.contourf(\n", " lon2D, lat2D, data2D, contours,\n", " cmap=cmap, vmin=vmin, vmax=vmax,\n", " transform=ccrs.PlateCarree())\n", "\n", " ax.add_feature(\n", " cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)\n", " ax.add_feature(cfeature.STATES, linewidth=0.5)\n", " ax.add_feature(cfeature.BORDERS, linewidth=0.5)\n", " ax.set_title(title, fontsize=fontSize)\n", "\n", " if cbar:\n", " cb = fig.colorbar(cf, ax=ax, format=labelFormat, pad=cPad)\n", " cb.ax.tick_params(labelsize=cFontSize)\n", " cb.set_label(cLabel, size=cFontSize)\n", "\n", " if cTicks is not None:\n", " tick_locator = ticker.MaxNLocator(nbins=cTicks)\n", " cb.locator = tick_locator\n", " cb.update_ticks()\n", "\n", " if saveFile == '':\n", " plt.show()\n", " else:\n", " plt.savefig(saveFile, bbox_inches=\"tight\")" ] }, { "cell_type": "markdown", "metadata": { "id": "tMJLENQNj4Yn" }, "source": [ "To make sure we have consistent colorscales between all our figures, we'll specify the colormap along with the `vmin` and `vmax` values of the colorabar based on the minimum and maximum values of the original swath data. We also need to set the extent of the area we want to show in the image, which we can do by specifying `DLAT` and `DLON` to specify the distance, in degrees, to show from the image center." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "J9tClF0nkHj5" }, "outputs": [], "source": [ "# Plot the satellite swath\n", "CMAP = 'viridis'\n", "VMIN = swath_h89.min()\n", "VMAX = swath_h89.max()\n", "DLAT = 5 # Extent of lat from center (degrees)\n", "DLON = 10 # Extent of lon from center (degrees)" ] }, { "cell_type": "markdown", "metadata": { "id": "qDre-WQpbb_6" }, "source": [ "The mapping function we have defined also allows us to specify the projection of the map, which does not have to be the same projection as the data, thanks to `cartopy`. Here we have chosen to show the data on the Robinson projection, to somewhat reduce the projection stretching in higher latitudes than would be for the Mercator projection." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "tS7LOqTcb9Hm" }, "outputs": [], "source": [ "PROJECTION = 'Robinson'" ] }, { "cell_type": "markdown", "metadata": { "id": "gT4lA53D_K0n" }, "source": [ "With these set, we can make our figure showing the original passive microwave swath data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 519 }, "executionInfo": { "elapsed": 30220, "status": "ok", "timestamp": 1686772295412, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "4KZoz2uA_Q_Y", "outputId": "08fc514c-7d42-4a35-e4bc-897dc2fec341" }, "outputs": [], "source": [ "map2D(swath_h89, swath_longitude, swath_latitude,\n", " cmap=CMAP,\n", " dLat=DLAT, dLon=DLON,\n", " projection=PROJECTION,\n", " vmin=VMIN, vmax=VMAX)" ] }, { "cell_type": "markdown", "metadata": { "id": "W8wsYqZUAPH6" }, "source": [ "Looking at this image, the horizontal polarization near 89 GHz shows ice scattering in deep convection surrounding the storm center as it travels on land. It also shows cirrus as well as warm rain against the low emissivity ocean background scene." ] }, { "cell_type": "markdown", "metadata": { "id": "wBdEZ3SA7CiS" }, "source": [ "## Regridding Swath Data\n", "When regridding swath data onto a uniform grid, you have to specify a projection. Multiple projections are available. Here, we explore regridding our swath data using the Mercator projection." ] }, { "cell_type": "markdown", "metadata": { "id": "6SwJ86vC7cVx" }, "source": [ "### Defining the swath\n", "In order for `pyresample` to properly interpret the satellite swath data, we have to define the area of the swath data. We can easily define this area using the swath latitude and longitude values from the TC PRIMED file.\n", "\n", "Recall that for `pyresample` we need to have the longitudes from -180 to 180, rather than 0 to 360 like that available in TC PRIMED files. While we have changed the longitude values when we loaded the data above, let's go ahead and check to make sure that our latitudes and longitudes are in the correct format for `pyresample` to interpret." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "arxxKPdt7Q-B" }, "outputs": [], "source": [ "# Latitude and longitude format check\n", "swath_longitude, swath_latitude = check_and_wrap(swath_longitude, swath_latitude)" ] }, { "cell_type": "markdown", "metadata": { "id": "cEfKHpg59qux" }, "source": [ "Now that we are sure that our longitude and latitude values are correctly formatted, we can define our satellite swath." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Eqo4-1AI7a9F" }, "outputs": [], "source": [ "# Swath definition\n", "swath = SwathDefinition(swath_longitude, swath_latitude)" ] }, { "cell_type": "markdown", "metadata": { "id": "msjP2b4aBFS1" }, "source": [ "### Setting up the Mercator projection\n", "In order to regrid the swath data onto a uniform grid, we first need to define the projection for the new grid. In `pyresample`, we do this by defining a dictionary where:\n", "* the projection is specificed as Mercator,\n", "* the central longitude is the longitude of the storm center\n", "* the radius for the projection is Earth's radius.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "06uF0CqWN-4X" }, "outputs": [], "source": [ "# Create a dictionary specifying the projection\n", "RADIUS_EARTH = 6371228.\n", "proj_dict_merc = {\n", " 'proj': 'merc',\n", " 'lon_0': storm_longitude,\n", " 'R': int(RADIUS_EARTH),\n", " 'units': 'm'}" ] }, { "cell_type": "markdown", "metadata": { "id": "9Gu0Rpu7Pkat" }, "source": [ "Next we can setup the specific area-related parameters for our Mercator projection.\n", "* To have the tropical cyclone be perfectly centered in our image, we'll create a new square grid, with matching `x` and `y` dimensions (*nx*, *ny*). Here, we'll create a new grid of 390 by 390 pixels.\n", "* Since we'd like to match other satellite data that may have higher resolutions, we'll use spatial resolutions (*dx*, *dy*) that correspond to 2 km. Using the great circle approximation, this evaluates to `dx=0.03593244699965226` and `dy=0.03593244699963805`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "LpliqM0KBBU3" }, "outputs": [], "source": [ "# Define the parameters necessary for the Mercator projection\n", "NX = 390 # number of pixels in x-direction (longitude)\n", "NY = 390 # number of pixels in y-direction (latitude)\n", "DX = 0.03593244699965226 # longitude spatial resolution (degrees)\n", "DY = 0.03593244699963805 # latitude spatial resolution (degrees)" ] }, { "cell_type": "markdown", "metadata": { "id": "xw3MDtNeFteP" }, "source": [ "Now that we have those parameters set, we can create our new projection area using `pyresample`'s `create_area_def` function. To keep the storm centered in our new projection, we will center the new projection on the storm longitude and latitude obtained from the TC PRIMED file." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "ES8RCy3xFsfn" }, "outputs": [], "source": [ "# Define the mercator projection\n", "AREA_ID = 'Mercator'\n", "merc_area = create_area_def(\n", " AREA_ID, proj_dict_merc,\n", " center=(storm_longitude, storm_latitude),\n", " shape=(NX, NY),\n", " resolution=(DX, DY),\n", " units='degrees')" ] }, { "cell_type": "markdown", "metadata": { "id": "_z5J8_2sIFSq" }, "source": [ "In order to plot the data on this new projection, we will need the associated latitude and longitude of this new grid area. Using `pyresample`'s area definition, we can easily pull out the full spatial 2D arrays.\n", "\n", "Since the Mercator projection is a uniform grid, we could reduce the arrays to 1D in order to minimize duplicity; however, for this notebook we will keep it in 2D to be consistent with the rest of TC PRIMED data. Note that if you do reduce down to 1D arrays, 2D data can be regenerated using [numpy's meshgrid function](https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "executionInfo": { "elapsed": 99, "status": "ok", "timestamp": 1686772912664, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "a9ShUZrHIsrB", "outputId": "14ff7c36-1f39-4a06-c895-a513d9a836f1" }, "outputs": [], "source": [ "# Pull out the latitude and longitude\n", "merc_longitude, merc_latitude = merc_area.get_lonlats()\n", "print(\"Mercator Longitude Shape: {}\".format(merc_longitude.shape))\n", "print(\"Mercator Latitude Shape: {}\".format(merc_latitude.shape))" ] }, { "cell_type": "markdown", "metadata": { "id": "F1jLQd1cRsyC" }, "source": [ "### How should you interpolate the data?\n", "To put the swath data onto the new grid, the data has to be interpolated to create new values corresponding to the new spatial locations. `pyresample` offers several different interpolation methods. Here, we demonstrate two [commonly used methods](https://gisgeography.com/raster-resampling/):\n", "* Nearest Neighbor\n", "* Bilinear Interpolation\n", "\n", "Both of these methods have different strengths and weaknesses. While we won't go into these methods in depth, here we're providing examples of both so that you can easily experiment with them to suit your needs." ] }, { "cell_type": "markdown", "metadata": { "id": "cgYfbZuoYPFd" }, "source": [ "#### Mercator projection using nearest neighbor\n", "Here, we'll use `pyresample` to regrid the swath data onto a Mercator grid using the nearest neighbor interpolation method. One important parameter of this method is the radius of influence, which is the radius around each grid pixel in meters to search for neighbours in the swath. Note that this distance needs to be large enough to ensure finding a neighbor for the target from the original data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "executionInfo": { "elapsed": 321, "status": "ok", "timestamp": 1686774058919, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "F8PKuQSFR_bH", "outputId": "610ecacb-320c-4e6f-d97c-7b597a9e9bd0" }, "outputs": [], "source": [ "# Nearest neighbor grid interpolation\n", "RADIUS_NN = 20000. # radius of influence (m)\n", "mercNN_h89 = kd_tree.resample_nearest(\n", " swath, swath_h89, merc_area,\n", " radius_of_influence=RADIUS_NN)\n", "print(\"Mercator Data Shape: {}\".format(mercNN_h89.shape))\n", "print(\" Min= {:.2f}, Max= {:.2f}\".format(\\\n", " mercNN_h89.min(), mercNN_h89.max()))" ] }, { "cell_type": "markdown", "metadata": { "id": "tN4owgWiZlrk" }, "source": [ "Printing the shape, we see as expected the new data has 390 by 390 pixels.\n", "\n", "`pyresample` provides a routine to quickly plot the data, which we can use as a way to quickly ensure there are no obvious errors in the new data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 406 }, "executionInfo": { "elapsed": 800, "status": "ok", "timestamp": 1686774072201, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "rQ0yD6X2ZocS", "outputId": "71f5f50e-a3d7-4b9c-b51d-0dd41b895089" }, "outputs": [], "source": [ "# Quickly plot the new data\n", "show_quicklook(merc_area, mercNN_h89,\n", " cmap=CMAP,\n", " vmin=VMIN, vmax=VMAX)" ] }, { "cell_type": "markdown", "metadata": { "id": "YNrBF-XAaqBB" }, "source": [ "To visually compare our results with the original swath data, we can use the same plotting function defined above to plot the data, using the latitude and longitude values from the Mercator area we defined." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 519 }, "executionInfo": { "elapsed": 46129, "status": "ok", "timestamp": 1686774126778, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "yXSjHPS7atXy", "outputId": "ca5c3f7d-afec-4199-9af0-c02d37c6a607" }, "outputs": [], "source": [ "# Plot the data on Mercator projection from 2D lat/lon\n", "map2D(mercNN_h89, merc_longitude, merc_latitude,\n", " cmap=CMAP,\n", " dLat=DLAT, dLon=DLON,\n", " projection=PROJECTION,\n", " vmin=VMIN, vmax=VMAX)" ] }, { "cell_type": "markdown", "metadata": { "id": "OrWvSY_Ac3ZY" }, "source": [ "Through visual comparison, the interpolated image look nearly identical to the swath image above!" ] }, { "cell_type": "markdown", "metadata": { "id": "AoyyKHp_c6tf" }, "source": [ "#### Mercator projection using bilinear interpolation\n", "Bilinear interpolation is another common interpolation method offered by `pyresample`. Key parameters for this method include:\n", "* radius: The radius around each target pixel in meters to search for neighbours in the source data. Again note that this distance needs to be large enough to ensure finding a neighbor for the target from the original data.\n", "* neighbours: The number of closest locations to consider when selecting the four data points around the target location. Note that this number needs to be large enough to ensure surrounding the target.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "CUeISHbDc4C1" }, "outputs": [], "source": [ "# Set bilinear interpolation parameters\n", "RADIUS_BILINEAR = 30000\n", "NEIGHBORS_BILINEAR = 32" ] }, { "cell_type": "markdown", "metadata": { "id": "R-lKmzdkd9xK" }, "source": [ "Similar to the steps above, we can convert the swath data to a Mercator projection using this interpolation method. This method requires you to first define a resampler object, where you specify the original area (here the swath area), the new area (here the Mercator grid), and the desired parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "FMNQPe3Xd82c" }, "outputs": [], "source": [ "bt_resampler = NumpyBilinearResampler(\n", " swath, merc_area,\n", " RADIUS_BILINEAR,\n", " neighbours=NEIGHBORS_BILINEAR)" ] }, { "cell_type": "markdown", "metadata": { "id": "jJbHf8NEf2zw" }, "source": [ "Once we have defined the resampler object, we can then use it on our swath data to perform the regridding." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "executionInfo": { "elapsed": 2258, "status": "ok", "timestamp": 1686774129426, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "oeVs4It9f_ZR", "outputId": "1bfc31c6-6116-4f65-9865-3d97bee011df" }, "outputs": [], "source": [ "mercB_h89 = bt_resampler.resample(swath_h89)\n", "print(\"Mercator Data Shape: {}\".format(mercB_h89.shape))\n", "print(\" Min= {:.2f}, Max= {:.2f}\".format(\\\n", " mercB_h89.min(), mercB_h89.max()))" ] }, { "cell_type": "markdown", "metadata": { "id": "FHOsLNCOh0u2" }, "source": [ "And now we can plot the results:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 519 }, "executionInfo": { "elapsed": 9206, "status": "ok", "timestamp": 1686774138630, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "uH4Vxs28gQjG", "outputId": "fc0e4c04-bffc-4aa6-c7ec-f56302631286" }, "outputs": [], "source": [ "# Plot the data on Mercator projection from 2D lat/lon\n", "map2D(mercB_h89, merc_longitude, merc_latitude,\n", " cmap=CMAP,\n", " dLat=DLAT, dLon=DLON,\n", " projection=PROJECTION,\n", " vmin=VMIN, vmax=VMAX)" ] }, { "cell_type": "markdown", "metadata": { "id": "om7cLR_EiKEe" }, "source": [ "Again we see results very similar to the original swath data, as well as from using nearest neighbor interpolation." ] }, { "cell_type": "markdown", "metadata": { "id": "5lljW5fwiNlA" }, "source": [ "### Lambert Azimuthal Equal Area (LAEA) projection\n", "As a final demonstration, we'll show how easily we can use `pyresample` to transform swath data to a different projection. Here we'll try the lambert azimuthal equal area projection, again centered on the tropical cyclone. Recall, the steps to do this are simple:\n", "1. Create swath definition (done!)\n", "2. Create new projection definition\n", "3. Specify the projection area\n", "4. Create new data using selected interpolation method" ] }, { "cell_type": "markdown", "metadata": { "id": "IDBRmP69jQT2" }, "source": [ "In defining the LAEA projection, we need to specify both the central latitude and the central longitude. Again, since we already have the storm location from TC PRIMED we will use these." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "JpYiU1fwh7e7" }, "outputs": [], "source": [ "# Define lambert azimuthal equal area projection\n", "laea_projection = {\n", " 'proj': 'laea',\n", " 'lat_0': storm_latitude,\n", " 'lon_0': storm_longitude,\n", " 'a': RADIUS_EARTH,\n", " 'units': 'm'}\n" ] }, { "cell_type": "markdown", "metadata": { "id": "WDkPzxkmjlJ_" }, "source": [ "In setting up the specific area for our LAEA projection, we will use parameter values to give us a similar area to our Mercator projection area.\n", "* To have the tropical cyclone be perfectly centered, we'll create a new square grid, with matching width and height dimensions with the same number of pixels as previously (*nx*, *ny*).\n", "* Rather than specifying the resolutions, this projection requires us to specify the area extent. Because we want a square with 2 km resolution, we set that up going equal distances from the center in both x- and y- directions." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "REbgJ0FIjnkh" }, "outputs": [], "source": [ "LAREA_WIDTH = NX # width size (pixels)\n", "LAREA_HEIGHT = NY # height size (pixels)\n", "LAREA_RES = 2000 # grid resolution (m)\n", "laea_extentx = LAREA_RES * LAREA_WIDTH # one-sided width extent\n", "laea_extenty = LAREA_RES * LAREA_HEIGHT # one-sided height extent\n", "laea_extent_area = (\\\n", " -laea_extentx, -laea_extenty,\n", " laea_extentx, laea_extenty)\n" ] }, { "cell_type": "markdown", "metadata": { "id": "Yes0UT-UmIEA" }, "source": [ "After setting the parameters, we can define the LAEA area using `pyresample`'s `AreaDefinition`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "XjZ538cBleLB" }, "outputs": [], "source": [ "LAREA_ID = 'LAEA'\n", "LAREA_DESCRIPTION = 'LAEA Example'\n", "LAREA_PROJECT_ID = 'TC-PRIMED Example'\n", "laea_area = AreaDefinition(\n", " LAREA_ID,\n", " LAREA_DESCRIPTION,\n", " LAREA_PROJECT_ID,\n", " laea_projection,\n", " LAREA_WIDTH,\n", " LAREA_HEIGHT,\n", " laea_extent_area)" ] }, { "cell_type": "markdown", "metadata": { "id": "t-4j1-XOmefi" }, "source": [ "Now we can interpolate the swath data, here using the nearest neighbors approach." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "executionInfo": { "elapsed": 424, "status": "ok", "timestamp": 1686774172705, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "m5xvdGA8mvhM", "outputId": "7f6fe24a-8268-4c43-821f-e66f68b31ade" }, "outputs": [], "source": [ "# Create regridded data using nearest neighbor\n", "laeaNN_h89 = kd_tree.resample_nearest(\n", " swath, swath_h89, laea_area,\n", " radius_of_influence=RADIUS_NN,\n", " fill_value=np.nan)\n", "print(\"LAEA Data Shape: {}\".format(laeaNN_h89.shape))\n", "print(\" Min= {:.2f}, Max= {:.2f}\".format(\\\n", " np.nanmin(laeaNN_h89), np.nanmax(laeaNN_h89)))" ] }, { "cell_type": "markdown", "metadata": { "id": "IkHV4cPTnXI_" }, "source": [ "Finally we can plot the results using latitudes and longitudes from the lambert equal area definition." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 519 }, "executionInfo": { "elapsed": 30350, "status": "ok", "timestamp": 1686774216289, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "YlvGCvzFnaAR", "outputId": "58142294-3ac6-4373-e2aa-844909ccc100" }, "outputs": [], "source": [ "# Plot lambert azimuthal equal area data\n", "laea_longitude, laea_latitude = laea_area.get_lonlats()\n", "map2D(laeaNN_h89, laea_longitude, laea_latitude,\n", " cmap=CMAP,\n", " dLat=DLAT, dLon=DLON,\n", " projection=PROJECTION,\n", " vmin=VMIN, vmax=VMAX)" ] }, { "cell_type": "markdown", "metadata": { "id": "ZfUNQfkjoe9v" }, "source": [ "Here we see that the values are similar, but the corners of the image are cutoff due to the different projection area." ] }, { "cell_type": "markdown", "metadata": { "id": "n_q9fPvjt3SR" }, "source": [ "## Comparing passive microwave data with precipitation\n", "Observations in the different passive microwave frequencies provide us with information on the various properties of the atmosphere (e.g., temperature, moisture, presence of liquid or frozen precipitation). These different information allow us to estimate the precipitation structure within the passive microwave swath. Along with passive microwave brightness temperature observations, TC PRIMED also contains passive-microwave-based rainfall retrievals from NASA's Goddard PROFiling algorithm (GPROF). Using the precipitation information from the same file as above, let's load, regrid and plot the GPROF surface precipitation rates from Hurricane Florence, and compare the surface precipitation rates with the passive microwave observation from above." ] }, { "cell_type": "markdown", "metadata": { "id": "LP_dRxE4T_Di" }, "source": [ "### Read and visualize precipitation data\n", "First, let's obtain the GPROF surface precipitation rate data from the same TC PRIMED overpass file as above. The surface precipitation rate data is stored in the `S1` swath of the GPROF group in the TC PRIMED overpass file." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "executionInfo": { "elapsed": 94, "status": "ok", "timestamp": 1686777637554, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "SspAWX56UDpY", "outputId": "df449030-5239-4220-e701-54f61872f52b" }, "outputs": [], "source": [ "# Read precipitation from the same TC PRIMED overpass file\n", "DS = Dataset(FILE_NAME, memory=url_response.content)\n", "precip_latitude = DS[\"GPROF/S1/latitude\"][:]\n", "precip_longitude = DS[\"GPROF/S1/longitude\"][:]\n", "precip_longitude[precip_longitude > 180] -= 360\n", "precip_data = DS[\"GPROF/S1/surfacePrecipitation\"][:]\n", "precip_units = DS[\"GPROF/S1/surfacePrecipitation\"].units\n", "print(\"GPROF Precipition Shape: {}\".format(precip_data.shape))\n", "print(\" Min= {:.2f}, Max= {:.2f} ({})\".\\\n", " format(precip_data.min(), precip_data.max(),\n", " precip_units))\n", "DS.close()" ] }, { "cell_type": "markdown", "metadata": { "id": "eZWJU-7dXc-X" }, "source": [ "To get an idea of what the precipitation looks like for Hurricane Florence at the overpass time, we can use our mapping routine to plot the precipitation data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 519 }, "executionInfo": { "elapsed": 4593, "status": "ok", "timestamp": 1686779453129, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "45cg6U4_VcZk", "outputId": "20bcd123-afb5-4e36-fc06-12571e9912b1" }, "outputs": [], "source": [ "# Map the original precipitation data\n", "map2D(precip_data, precip_longitude, precip_latitude,\n", " cLabel='mm/h', cmap=CMAP,\n", " dLat=DLAT, dLon=DLON,\n", " projection=PROJECTION)" ] }, { "cell_type": "markdown", "metadata": { "id": "Eoyk8Y9uVzwt" }, "source": [ "Looking closely at the map of the original precipitation swath, notice that there are points over western Pennsylvania and southwestern New York where data is missing. We can use the data mask to set any missing value to zero to ensure that the regridding process goes smoothly." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "NOy40P1H7wfw" }, "outputs": [], "source": [ "# Replace any bad values with zero\n", "precip_data = np.where(precip_data.mask, 0.0, precip_data)" ] }, { "cell_type": "markdown", "metadata": { "id": "SfjQ64zbbPQr" }, "source": [ "### Regrid the precipitation data\n", "Next, in order to directly compare the passive microwave data with the precipitation data, we can put them both on the same uniform grid. Here we will use the Mercator grid that we have already used for the microwave data. To regrid the precipitation data, recall we first need to define the swath area of the precipitation data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "9uLMTZQpXCGw" }, "outputs": [], "source": [ "# Define precipitation swath area\n", "precip_swath = SwathDefinition(precip_longitude, precip_latitude)" ] }, { "cell_type": "markdown", "metadata": { "id": "NAb6NhivWHYu" }, "source": [ "Using the swath areal definition, now we can interpolate the precipitation swath onto the Mercator grid." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "executionInfo": { "elapsed": 274, "status": "ok", "timestamp": 1686779890857, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "MWYGy6kiXxke", "outputId": "032cfa8d-6258-4b17-e2cf-ba946c6d7231" }, "outputs": [], "source": [ "# Regrid precipitation onto Mercator grid\n", "# using nearest neighbor interpolation\n", "mercNN_precip = kd_tree.resample_nearest(\n", " precip_swath, precip_data, merc_area,\n", " fill_value=0.,\n", " radius_of_influence=RADIUS_NN)\n", "print(\"Regridded precip Min={:.2f}, Max={:.2f}\".\\\n", " format(mercNN_precip.min(), mercNN_precip.max()))" ] }, { "cell_type": "markdown", "metadata": { "id": "e4b8MRS2mbTJ" }, "source": [ "Just for a sanity check, let's quickly check the size of our regridded precipitation data with the regridded 89-GHz passive microwave data. Since they were both regridded to the same Mercator area that we defined (merc_area), they should both have the same shape." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "executionInfo": { "elapsed": 105, "status": "ok", "timestamp": 1686779908611, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "jgu08-K0U9ly", "outputId": "5b0e1715-abe8-48bb-97f2-37a360ce715c" }, "outputs": [], "source": [ "print(f\"Size of 89-GHz regridded: {mercNN_h89.shape}\")\n", "print(f\"Size of precip regridded: {mercNN_precip.shape}\")" ] }, { "cell_type": "markdown", "metadata": { "id": "QNM3UjPsm_VN" }, "source": [ "### Compare the passive microwave and precipitation data\n", "\n", "Now that we have the passive microwave and precipitation data on the same grid, we can plot them next to each other to see what their images look like for Hurricane Florence. To do this, let's create a plotting routine to show two images next to each other." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Zq4Ocg-5fJ05" }, "outputs": [], "source": [ "def imageArea2(areaDef, data1, data2,\n", " cbar=True,\n", " cbarLabel1='K',\n", " cbarLabel2='mm/h',\n", " cbarFontSize=10,\n", " cbarPad = 0.01,\n", " cbarShrink=0.45,\n", " cmap1='viridis',\n", " cmap2='viridis',\n", " origin='upper',\n", " title1='',\n", " title2='',\n", " titleFontSize=12,\n", " vmin1=None,\n", " vmax1=None,\n", " vmin2=None,\n", " vmax2=None\n", " ):\n", "\n", " \"\"\"\n", " Create a figure with two side-by-side images from area definitions.\n", "\n", " Parameters\n", " ----------\n", " areaDef : area definition\n", " Pyresample AreaDefinition that can be used to create\n", " a corresponding cartopy coordinate reference\n", " data1 : array-like\n", " array of 2D data matching the area definition,\n", " which will be plotted in the left-hand image\n", " data2 : array-like\n", " array of 2D data matching the area definition,\n", " which will be plotted in the right-hand image\n", " cbar : bool, default=True\n", " flag to show colorbars\n", " cbarLabel1 : str, default='K'\n", " colorbar label for data1\n", " cbarLabel2 : str, default='mm/h'\n", " colorbar label for data2\n", " cbarFontSize : int, default=10\n", " colorbar label font size\n", " cbarPad : float, default= 0.01\n", " colorbar padding\n", " cbarShrink : float, default=0.45\n", " shrink ratio for colorbar\n", " cmap1 : str, default='viridis'\n", " color table for data1\n", " cmap2 : str, default='viridis'\n", " color table for data2\n", " origin : str, default='upper'\n", " orientation of area\n", " title1 : str, default=''\n", " figure title for data1\n", " title2 : str, default=''\n", " title for data2\n", " titleFontSize : int, default=12\n", " font size for titles\n", " vmin1 : float, optional\n", " minimum value for data1 colorbar, calculated from data if None\n", " vmax1 : float, optional\n", " maximum value for data1 colorbar, calculated from data if None\n", " vmin2 : float, optional\n", " minimum value for data2 colorbar, calculated from data if None\n", " vmax2 : float, optional\n", " maximum value for data2 colorbar, calculated from data if None\n", "\n", " Yields\n", " ------\n", " side-by-side image of two regridded satellite swath data on a map\n", " with a projection\n", " \"\"\"\n", "\n", " fig = plt.figure\n", "\n", " if vmin1 is None:\n", " vmin1 = np.nanmin(data1)\n", " if vmax1 is None:\n", " vmax1 = np.nanmax(data1)\n", " if vmin2 is None:\n", " vmin2 = np.nanmin(data2)\n", " if vmax2 is None:\n", " vmax2 = np.nanmax(data2)\n", "\n", " # Create the figure and plot the background\n", " crs = areaDef.to_cartopy_crs()\n", " extent = crs.bounds\n", " fig, ax = plt.subplots(\\\n", " ncols=2, nrows=1,\n", " subplot_kw={'projection': crs})\n", "\n", " cbarLList = [cbarLabel1, cbarLabel2]\n", " cmapList = [cmap1, cmap2]\n", " dataList = [data1, data2]\n", " titleList = [title1, title2]\n", " vminList = [vmin1, vmin2]\n", " vmaxList = [vmax1, vmax2]\n", " for i in range(2):\n", " im = ax[i].imshow(\\\n", " dataList[i],\n", " extent=extent,\n", " origin=origin,\n", " cmap=cmapList[i],\n", " vmin=vminList[i],\n", " vmax=vmaxList[i],\n", " transform=crs)\n", " ax[i].coastlines(resolution='50m',\n", " color='black',\n", " linewidth=0.5)\n", " ax[i].add_feature(cfeature.STATES, linewidth=0.5)\n", " ax[i].add_feature(cfeature.BORDERS, linewidth=0.5)\n", " ax[i].set_title(titleList[i], fontsize=titleFontSize)\n", "\n", " if cbar:\n", " cb = fig.colorbar(\\\n", " im, ax=ax[i],\n", " pad=cbarPad,\n", " shrink=cbarShrink)\n", " cb.set_label(cbarLList[i], size=cbarFontSize)\n", " cb.ax.tick_params(labelsize=cbarFontSize)\n", "\n", " plt.tight_layout()\n", " plt.show()\n", " plt.close()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 335 }, "executionInfo": { "elapsed": 14684, "status": "ok", "timestamp": 1686779976176, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "urIOJO9-jtur", "outputId": "6add378c-65da-45b2-b1f4-e12941ff8611" }, "outputs": [], "source": [ "imageArea2(merc_area, mercNN_h89, mercNN_precip,\n", " title1='89 GHz', title2='Precipitation')" ] }, { "cell_type": "markdown", "metadata": { "id": "ENjnQOh9lMUS" }, "source": [ "In the 89 GHz passive microwave data, we can see the eye of Hurricane Florence with high brightness temperatures, surrounded by spiralling areas of low brightness temperatures associated with deep convection. Further from the storm, we also see bands of low brightness temperatures just off the coasts of North Carolina and Virginia in the Atlantic. In the corresponding precipitation data, we see intense precipitation in the south-eastern quadrant of the storm, as well as heavy precipitation that appears to coincide with the low 89 GHz brightness temperatures off the coast in the Atlantic." ] }, { "cell_type": "markdown", "metadata": { "id": "IZxaysl3mpun" }, "source": [ "Because the datasets are now on the same grid, we can actually overlay the precipitation data on the passive microwave data to see how the two line up. Let's again define a routine and plot this up." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "VEdEKrett4-z" }, "outputs": [], "source": [ "def imageArea_overlay(\n", " areaDef,\n", " data1, data2,\n", " cbar1=True,\n", " cbar2=True,\n", " cbarAspect=50,\n", " cbarFontSize=10,\n", " cbarLabel1='K',\n", " cbarLabel2='mm/h',\n", " cbarPad1=0.03,\n", " cbarPad2=0.01,\n", " cbarShrink1=0.65,\n", " cbarShrink2=0.45,\n", " cmap1='viridis',\n", " cmap2='Greys',\n", " coast_res='50m',\n", " figSize=(6, 6),\n", " marker='x',\n", " npts=80,\n", " nlevs=101,\n", " origin='upper',\n", " title='',\n", " titleSize=14,\n", " vmin1=None,\n", " vmax1=None,\n", " vmin2=None,\n", " vmax2=None):\n", "\n", " \"\"\"\n", " Create a figure with an image of one dataset and\n", " overlaying symbols colored to match a second dataset.\n", "\n", " Parameters\n", " ----------\n", " areaDef : area definition\n", " Pyresample AreaDefinition that can be used to create\n", " a corresponding cartopy coordinate reference\n", " data1 : array-like\n", " array of 2D data matching the area definition,\n", " which will be plotted in the left-hand image\n", " data2 : array-like\n", " array of 2D data matchin the area definition,\n", " which will be plotted in the right-hand image\n", " cbar1 : bool, default=True\n", " flag to show colorbar for data1\n", " cbar2 : bool, default=True\n", " flag to show colorbar for data2\n", " cbarAspect : int, default=50\n", " colorbar aspect ratio\n", " cbarFontSize : int, default=10\n", " colorbar font size\n", " cbarLabel1 : str, default='K'\n", " colorbar label for data1\n", " cbarLabel2 : str. default='mm/h'\n", " colorbar label for data2\n", " cbarPad1 : float, default=0.03\n", " padding for data1 colorbar\n", " cbarPad2 : float, default=0.01\n", " padding for data2 colorbar\n", " cbarShrink1 : float, default=0.65\n", " shrink ratio for data1 colorbar\n", " cbarShrink2 : float, default=0.45\n", " shrink ratio for data2 colorbar\n", " cmap1 : str, default='viridis'\n", " color table for data1\n", " cmap2 : str, default='Greys'\n", " color table for data2\n", " coast_res : str, default='50m'\n", " resolution of coast overlay\n", " figSize : tuple, default=(6, 6)\n", " figure size\n", " marker : str, default='x'\n", " marker style for the data2 overlay\n", " npts : int, default=80\n", " number of points of data2 to overlay\n", " nlevs : int, default=101\n", " number of levels in contours\n", " origin : str, default='upper'\n", " orientation of the data images\n", " title : str, default=''\n", " title for the image\n", " titleSize : int, default=14\n", " font size of the title\n", " vmin1 : float, optional\n", " minimum value for data1 colorbar, calculated from data if None\n", " vmax1 : float, optional\n", " maximum value for data1 colorbar, calculated from data if None\n", " vmin2 : float, optional\n", " minimum value for data2 colorbar, calculated from data if None\n", " vmax2 : float, optional\n", " maximum value for data2 colorbar, calculated from data if None\n", "\n", " Yields\n", " ------\n", " image of regridded satellite swath data on a map with a projection\n", " with data2 overlaid on data1\n", " \"\"\"\n", "\n", " if vmin1 is None:\n", " vmin1 = np.nanmin(data1)\n", " if vmax1 is None:\n", " vmax1 = np.nanmax(data1)\n", " if vmin2 is None:\n", " vmin2 = np.nanmin(data2)\n", " if vmax2 is None:\n", " vmax2 = np.nanmax(data2)\n", "\n", " # setup coordinates\n", " crs = areaDef.to_cartopy_crs()\n", " extent = crs.bounds\n", "\n", " # create figure\n", " fig = plt.figure(figsize=figSize)\n", " ax = plt.axes(projection=crs)\n", " ax.set_extent(extent, crs)\n", "\n", " # add first image\n", " cf = ax.imshow(data1, transform=crs,\n", " extent=extent,\n", " origin=origin,\n", " cmap=cmap1, vmin=vmin1, vmax=vmax1)\n", " ax.coastlines(resolution=coast_res, color='black', linewidth=0.5)\n", " ax.add_feature(cfeature.STATES, linewidth=0.5)\n", " ax.add_feature(cfeature.BORDERS, linewidth=0.5)\n", " if title:\n", " ax.set_title(title, fontsize=titleSize)\n", "\n", " # add overlay\n", " nXY = data1.shape[-1]\n", " ptsArray = np.linspace(0, nXY - 1, npts, dtype=np.int32)\n", "\n", " aArray = np.ones((nlevs))\n", " if nlevs > 40:\n", " aArray[:10] = 0.1\n", " aArray[10:20] = 0.3\n", " aArray[20:30] = 0.5\n", " aArray[30:40] = 0.7\n", " aArray[40:] = 0.8\n", " sizeArray = np.linspace(3, 4, nlevs)\n", " xArray = np.linspace(0.0, 1.0, nlevs)\n", "\n", " lons, lats = areaDef.get_lonlats()\n", " range2 = vmax2 - vmin2\n", " rgb = cm.get_cmap(cmap2)(xArray)[:, :3]\n", " for j in ptsArray:\n", " for i in ptsArray:\n", " subHere = np.min(\\\n", " [np.max([data2[j, i] - vmin2, 0.]),\n", " range2])\n", " valHere = int(subHere / range2 * 100.)\n", " colHere = rgb[valHere, :]\n", " alphaHere = aArray[valHere]\n", " sizeHere = sizeArray[valHere]\n", " ax.plot(lons[j, i], lats[j, i],\n", " marker,\n", " alpha=alphaHere,\n", " color=colHere,\n", " markersize=sizeHere,\n", " transform=ccrs.Geodetic())\n", "\n", " if cbar1:\n", " cb = fig.colorbar(\\\n", " cf, ax=ax,\n", " pad=cbarPad1,\n", " shrink=cbarShrink1)\n", " cb.set_label(cbarLabel1, size=cbarFontSize)\n", "\n", " if cbar2:\n", " norm = Normalize(vmin=vmin2, vmax=vmax2)\n", " cbu = fig.colorbar(\n", " cm.ScalarMappable(norm=norm, cmap=cmap2),\n", " ax=ax, shrink=cbarShrink2, pad=cbarPad2)\n", " cbu.set_label(cbarLabel2, size=cbarFontSize)\n", "\n", " plt.tight_layout()\n", " plt.show()\n", " plt.close()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 510 }, "executionInfo": { "elapsed": 12644, "status": "ok", "timestamp": 1686779998299, "user": { "displayName": "Naufal Razin", "userId": "00945304052336708650" }, "user_tz": 360 }, "id": "VIgamGsMxwZY", "outputId": "60c677b5-67ab-4fcc-ab7f-b0fa9d59a5d5" }, "outputs": [], "source": [ "imageArea_overlay(merc_area, mercNN_h89, mercNN_precip,\n", " vmax2=20)" ] }, { "cell_type": "markdown", "metadata": { "id": "Ax0gkfTgoXjQ" }, "source": [ "Here we see that heavier precipitation does coincide with the low brightness temperatures for the 89 GHz microwave data. Particularly, we see intense precipitation to the south-east of the eye covering the spiral band of low brightness temperatures to the south of the eye, as well as intense precipitation occurring where there are low-brightness temperatures seen off the North Carolina and Virgina coasts." ] }, { "cell_type": "markdown", "metadata": { "id": "XwS5tFCwp9lj" }, "source": [ "## Final thoughts\n", "When applying machine learning methods to the passive microwave data from TC PRIMED, you may want to pre-process the swath data by regridding the swath onto a common, uniform grid. One way to regrid the data is by using existing Python packages like `pyresample`. In this notebook, we have demonstrated examples of how to use `pyresample` to regrid TC PRIMED passive microwave data onto a grid using two different types of projections (Mercator and Lambert Azimuthal Equidistant Area) and two different types of interpolation methods (nearest neighbor and bilinear). We then demonstrated how this could be useful to overlay different datasets, such as passive microwave data and precipitation. These examples serve as guidance for data pre-processing when applying more complex machine learning approaches, like convolutional neural networks, to TC PRIMED data." ] }, { "cell_type": "markdown", "metadata": { "id": "Elm502sjrYHE" }, "source": [ "## Practice\n", "1. Projection resolution\n", " * Change the resolution size for the projections (*i.e.*, for Mercator change dx/dy and for LAEA change lea_extent). How does that change the projected image?\n", "\n", "\n", "2. Projection size\n", " * Change the size of the projected data (*i.e.*, change nx/ny). How does this change the projection differently than changing the resolution?\n", "\n", "\n", "3. Data variable\n", " * Read in different swath data available in TC PRIMED (*i.e.*, different polarization or frequency). See how it correlates with additional sources of data, and think through how regridding would be useful for collocating data, such as for an ML application." ] }, { "cell_type": "markdown", "metadata": { "id": "w2BB4VMJqZT7" }, "source": [ "## Data Statement\n", "- Razin, Muhammad Naufal; Slocum, Christopher J.; Knaff, John A.; Brown, Paula J. 2023. Tropical Cyclone PRecipitation, Infrared, Microwave, and Environmental Dataset (TC PRIMED). v01r00. NOAA National Centers for Environmental Information. https://doi.org/10.25921/dmy1-0595." ] }, { "cell_type": "markdown", "metadata": { "id": "DNU21CXFtcz1" }, "source": [ "## References\n", "Material in this notebook is taken from [Pyresample documentation](https://pyresample.readthedocs.io/en/latest/index.html). Pyresample is available on GitHub at https://github.com/pytroll/pyresample.\n" ] }, { "cell_type": "markdown", "metadata": { "id": "RWOEHJcUqs2t" }, "source": [ "## Metadata\n", "* Language / package(s)\n", " * Python\n", " * NetCDF\n", " * Matplotlib\n", " * Cartopy\n", " * Pyresample\n", "* Domain\n", " * NOAA\n", " * NASA\n", " * JAXA\n", "* Application keywords\n", " * Satellite passive microwave\n", " * Resampling\n", " * Data collocation\n", "* Geophysical keywords\n", " * Satellite data\n", " * Tropical Cyclones\n", "* AI keywords\n", " * Data pre-processing" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "rILbAJ9aqZ-t" }, "outputs": [], "source": [] } ], "metadata": { "colab": { "provenance": [ { "file_id": "1HXefT9zZdoIBP-WlDW762EbM4eYJd9pF", "timestamp": 1686563716335 } ] }, "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.9.15" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 4 }