Vorticity_Divergence_Inversion.ipynb

{"cells": [{"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import os.path\n", "import xarray as xr\n", "import numpy as np\n", "import boto3\n", "import metpy.calc as mpcalc\n", "from botocore import UNSIGNED\n", "from botocore.config import Config\n", "import matplotlib.pyplot as plt\n", "from matplotlib.cm import get_cmap\n", "#from matplotlib.colormaps import get_cmap\n", "import matplotlib.ticker as mticker\n", "import cartopy.crs as crs\n", "from cartopy.feature import NaturalEarthFeature\n", "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n", "import sys\n", "import time \n", "if (not os.path.isfile('gfs.t12z.pgrb2.0p25.f000')):\n", "    client = boto3.client('s3', config=Config(signature_version=UNSIGNED))\n", "    client.download_file('noaa-gfs-bdp-pds', 'gfs.20230809/12/atmos/gfs.t12z.pgrb2.0p25.f000', 'gfs.t12z.pgrb2.0p25.f000')"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["u850 = xr.open_dataset('gfs.t12z.pgrb2.0p25.f000', engine='cfgrib',backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa', 'shortName': 'u', 'level': 850}})\n", "u = u850.u\n", "print(u.shape)\n", "v850 = xr.open_dataset('gfs.t12z.pgrb2.0p25.f000', engine='cfgrib', backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa', 'shortName': 'v', 'level': 850}})\n", "v = v850.v"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the 850 hPa relative vorticity."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["vort850 = mpcalc.vorticity(u, v)\n", "fig = plt.figure(figsize=(12,9), dpi=300.)\n", "# Create a set of axes for the figure and set\n", "# its map projection to that of the input data.\n", "ax = plt.axes(projection=crs.PlateCarree())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add country borders and coastlines."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["countries = NaturalEarthFeature(category=\"cultural\", scale=\"50m\",\n", "                                      facecolor=\"none\",\n", "                                      name=\"admin_0_countries\")\n", "ax.add_feature(countries, linewidth=.5, edgecolor=\"black\")\n", "ax.coastlines('50m', linewidth=0.8)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plot = vort850.plot(levels=np.arange(-1.e-4, 1.e-4, 0.2e-5), cmap=get_cmap('PRGn'), transform=crs.PlateCarree(), cbar_kwargs={'label':'relative vorticity (x$10^{-5} s^{-1}$)', 'shrink': 0.98})"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Set the map's extent to cover just Hurricane Dora."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["ax.set_extent([-180.,-150.,0.,20.],crs=crs.PlateCarree())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add latitude/longitude gridlines."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["gridlines = ax.gridlines(color=\"grey\", linestyle=\"dotted\", draw_labels=True)\n", "gridlines.xlabels_top = False\n", "gridlines.ylabels_right = False\n", "gridlines.xlocator = mticker.FixedLocator(np.arange(-180.,149.,5.))\n", "gridlines.ylocator = mticker.FixedLocator(np.arange(0.,21.,5.))\n", "gridlines.xlabel_style = {'size':12, 'color':'black'}\n", "gridlines.ylabel_style = {'size':12, 'color':'black'}\n", "gridlines.xformatter = LONGITUDE_FORMATTER\n", "gridlines.yformatter = LATITUDE_FORMATTER"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add a plot title, then show the image."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plt.title(\"GFS 0-h 850 hPa relative vorticity (x$10^{-5} s^{-1}$) at 1200 UTC 9 August 2023\")\n", "plt.show()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the 850 hPa divergence."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["div850 = mpcalc.divergence(u, v)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Create a figure instance."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["fig = plt.figure(figsize=(12,9), dpi=300.)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Create a set of axes for the figure and set<br>\n", "its map projection to that of the input data."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["ax = plt.axes(projection=crs.PlateCarree())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add country borders and coastlines."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["countries = NaturalEarthFeature(category=\"cultural\", scale=\"50m\",\n", "                                      facecolor=\"none\",\n", "                                      name=\"admin_0_countries\")\n", "ax.add_feature(countries, linewidth=.5, edgecolor=\"black\")\n", "ax.coastlines('50m', linewidth=0.8)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Plot the 850 hPa divergence using xarray's plot functionality."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plot = div850.plot(levels=np.arange(-1.e-4, 1.e-4, 0.2e-5), cmap=get_cmap('PRGn'), transform=crs.PlateCarree(), cbar_kwargs={'label':'relative vorticity (x$10^{-5} s^{-1}$)', 'shrink': 0.98})"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Set the map's extent to cover just Hurricane Dora."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["ax.set_extent([-180.,-150.,0.,20.],crs=crs.PlateCarree())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add latitude/longitude gridlines."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["gridlines = ax.gridlines(color=\"grey\", linestyle=\"dotted\", draw_labels=True)\n", "gridlines.xlabels_top = False\n", "gridlines.ylabels_right = False\n", "gridlines.xlocator = mticker.FixedLocator(np.arange(-180.,149.,5.))\n", "gridlines.ylocator = mticker.FixedLocator(np.arange(0.,21.,5.))\n", "gridlines.xlabel_style = {'size':12, 'color':'black'}\n", "gridlines.ylabel_style = {'size':12, 'color':'black'}\n", "gridlines.xformatter = LONGITUDE_FORMATTER\n", "gridlines.yformatter = LATITUDE_FORMATTER"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add a plot title, then show the image."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plt.title(\"GFS 0-h 850 hPa divergence (x$10^{-5} s^{-1}$) at 1200 UTC 9 August 2023\")\n", "plt.show()"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["vortmask = mpcalc.bounding_box_mask(vort850,5.,13.5,191.,202.)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["divmask = mpcalc.bounding_box_mask(div850,5.,13.5,191.,202.)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["i_bb_indices = mpcalc.find_bounding_box_indices(vortmask,5.,13.5,191.,202.)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["o_bb_indices = mpcalc.find_bounding_box_indices(vortmask,0.,30,180.,220)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["dx, dy = mpcalc.lat_lon_grid_deltas(vortmask.longitude, vortmask.latitude)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["upsi,vpsi  = mpcalc.rotational_wind_from_inversion(vortmask,dx,dy,o_bb_indices,i_bb_indices)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Create a figure instance."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["fig = plt.figure(figsize=(12,9), dpi=300.)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Create a set of axes for the figure and set<br>\n", "its map projection to that of the input data."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["ax = plt.axes(projection=crs.PlateCarree())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add country borders and coastlines."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["countries = NaturalEarthFeature(category=\"cultural\", scale=\"50m\",\n", "                                      facecolor=\"none\",\n", "                                      name=\"admin_0_countries\")\n", "ax.add_feature(countries, linewidth=.5, edgecolor=\"black\")\n", "ax.coastlines('50m', linewidth=0.8)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the magnitude of the non-divergent component of the 850 hPa wind."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["nd_spd = np.sqrt(upsi**2 + vpsi**2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Plot this using xarray's plot functionality."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plot = nd_spd.plot(levels=np.arange(0., 13., 1.), cmap=get_cmap('YlGnBu'), transform=crs.PlateCarree(), cbar_kwargs={'label':'non-divergent wind ($m s^{-1}$)', 'shrink': 0.98})"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Set the map's extent to match that over which we computed the non-divergent wind."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["ax.set_extent([-180.,-140.,0.,30.],crs=crs.PlateCarree())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add latitude/longitude gridlines."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["gridlines = ax.gridlines(color=\"grey\", linestyle=\"dotted\", draw_labels=True)\n", "gridlines.xlabels_top = False\n", "gridlines.ylabels_right = False\n", "gridlines.xlocator = mticker.FixedLocator(np.arange(-180.,139.,5.))\n", "gridlines.ylocator = mticker.FixedLocator(np.arange(0.,31.,5.))\n", "gridlines.xlabel_style = {'size':12, 'color':'black'}\n", "gridlines.ylabel_style = {'size':12, 'color':'black'}\n", "gridlines.xformatter = LONGITUDE_FORMATTER\n", "gridlines.yformatter = LATITUDE_FORMATTER"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add a plot title, then show the image."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plt.title(\"GFS 0-h 850 hPa non-divergent wind magnitude ($m s^{-1}$) due to Dora at 1200 UTC 9 August 2023\")\n", "plt.show()"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["uchi,vchi = mpcalc.divergent_wind_from_inversion(divmask,dx,dy,o_bb_indices,i_bb_indices)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Create a set of axes for the figure and set<br>\n", "its map projection to that of the input data."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["ax = plt.axes(projection=crs.PlateCarree())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add country borders and coastlines."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["countries = NaturalEarthFeature(category=\"cultural\", scale=\"50m\",\n", "                                      facecolor=\"none\",\n", "                                      name=\"admin_0_countries\")\n", "ax.add_feature(countries, linewidth=.5, edgecolor=\"black\")\n", "ax.coastlines('50m', linewidth=0.8)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the magnitude of the non-divergent component of the 850 hPa wind."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["nd_spd = np.sqrt(uchi**2 + vchi**2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Plot this using xarray's plot functionality."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plot = nd_spd.plot(levels=np.arange(0., 13., 1.), cmap=get_cmap('YlGnBu'), transform=crs.PlateCarree(), cbar_kwargs={'label':'non-divergent wind ($m s^{-1}$)', 'shrink': 0.98})"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Set the map's extent to match that over which we computed the non-divergent wind."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["ax.set_extent([-180.,-140.,0.,30.],crs=crs.PlateCarree())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add latitude/longitude gridlines."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["gridlines = ax.gridlines(color=\"grey\", linestyle=\"dotted\", draw_labels=True)\n", "gridlines.top_labels = False\n", "gridlines.right_labels = False\n", "gridlines.xlocator = mticker.FixedLocator(np.arange(-180.,139.,5.))\n", "gridlines.ylocator = mticker.FixedLocator(np.arange(0.,31.,5.))\n", "gridlines.xlabel_style = {'size':12, 'color':'black'}\n", "gridlines.ylabel_style = {'size':12, 'color':'black'}\n", "gridlines.xformatter = LONGITUDE_FORMATTER\n", "gridlines.yformatter = LATITUDE_FORMATTER"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Add a plot title, then show the image."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plt.title(\"GFS 0-h 850 hPa divergent wind magnitude ($m s^{-1}$) due to Dora at 1200 UTC 9 August 2023\")\n", "plt.savefig('vectorized_version')\n", "plt.show()"]}], "metadata": {"kernelspec": {"display_name": "Python 3", "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.6.4"}}, "nbformat": 4, "nbformat_minor": 2}