{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MetPy Intro: NYS Mesonet Map\n", "---" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Overview\n", "In this notebook, we'll use Cartopy, Matplotlib, and Pandas (with a bunch of help from [MetPy](https://unidata.github.io/MetPy)) to read in, manipulate, and visualize current data from the [New York State Mesonet](https://www2.nysmesonet.org)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prerequisites\n", "\n", "| Concepts | Importance | Notes |\n", "| --- | --- | --- |\n", "| Matplotlib | Necessary | |\n", "| Cartopy | Necessary | |\n", "| Pandas | Necessary | |\n", "| MetPy | Necessary | Intro |\n", "\n", "* **Time to learn**: 30 minutes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "___\n", "## Imports" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "from datetime import datetime\n", "from cartopy import crs as ccrs\n", "from cartopy import feature as cfeature\n", "from metpy.calc import wind_components, dewpoint_from_relative_humidity\n", "from metpy.units import units\n", "from metpy.plots import StationPlot, USCOUNTIES" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a Pandas `DataFrame` object pointing to the most recent set of NYSM obs." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nysm_data = pd.read_csv('https://www.atmos.albany.edu/products/nysm/nysm_latest.csv')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nysm_data.columns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create several `Series` objects for some of the columns." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stid = nysm_data['station']\n", "lats = nysm_data['lat']\n", "lons = nysm_data['lon']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our goal is to make a map of NYSM observations, which includes the wind velocity. The convention is to plot wind velocity using wind barbs. The **MetPy** library allows us to not only make such a map, but perform a variety of meteorologically-relevant calculations and diagnostics. Here, we will use such a calculation, which will determine the two scalar components of wind velocity (*u* and *v*), from wind speed and direction. We will use MetPy's [wind_components](https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.wind_components.html) method. \n", "\n", "This method requires us to do the following:\n", "1. Create Pandas `Series` objects for the variables of interest\n", "1. Extract the underlying Numpy array via the `Series`' `values` attribute\n", "1. Attach units to these arrays using MetPy's `units` class" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Perform these three steps" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "wspd = nysm_data['max_wind_speed_prop [m/s]'].values * units['m/s']\n", "drct = nysm_data['wind_direction_prop [degrees]'].values * units['degree']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Tip: The * operator might make you think that we're multiplying. However, in this case we use it to \"attach\" units to an array of values.
" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "Examine these two *units aware* `Series`" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "wspd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Tip: This is a unique type of object: not only does it have values, aka magnitude, but it also has Units attached. This is a very nice property, since it will make calculations that require awareness of units much easier!
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "drct" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convert wind speed from m/s to knots" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "wspk = wspd.to('knots')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "wspk" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Perform the vector decomposition" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u, v = wind_components(wspk, drct)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Take a look at one of the output components:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Exercise: Create a 2-meter temperature object; attach units of degrees Celsius (degC) to it, and then convert to Fahrenheit (degF).
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Write your code here" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [] }, "outputs": [], "source": [ "#%load /spare11/atm350/common/mar12/metpy1.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's plot several of the meteorological values on a map. We will use **Matplotlib** and **Cartopy**, as well as **MetPy**'s `StationPlot` method." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create units-aware objects for relative humidity and station pressure" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "rh = nysm_data['relative_humidity [percent]'].values * units('percent')\n", "pres = nysm_data['station_pressure [mbar]'].values * units('mbar')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the map, centered over NYS, with add some geographic features, and the mesonet data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Determine the current time, for use in the plot's title and a version to be saved to disk." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "timeString = nysm_data['time'][0]\n", "timeObj = datetime.strptime(timeString,\"%Y-%m-%d %H:%M:%S\")\n", "titleString = datetime.strftime(timeObj,\"%B %d %Y, %H%M UTC\")\n", "figString = datetime.strftime(timeObj,\"%Y%m%d_%H%M\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Previously, we used the `ax.scatter` and `ax.text` methods to plot markers for the stations and their site id's. The latter method does not provide an intuitive means to orient several text stings relative to a central point as we typically do for a meteorological surface station plot. Let's take advantage of the `Metpy` package, and its StationPlot method!\n", "\n", "#### Be patient: this may take a minute or so to plot, if you chose the highest resolution for the shapefiles!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Specify two resolutions: one for the Natural Earth shapefiles, and the other for MetPy's US County shapefiles." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = '10m'\n", "resCounty = '5m'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Set the domain for defining the plot region.\n", "latN = 45.2\n", "latS = 40.2\n", "lonW = -80.0\n", "lonE = -72.0\n", "cLat = (latN + latS)/2\n", "cLon = (lonW + lonE )/2\n", "\n", "proj = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)\n", "\n", "fig = plt.figure(figsize=(18,12),dpi=150) # Increase the dots per inch from default 100 to make plot easier to read\n", "ax = fig.add_subplot(1,1,1,projection=proj)\n", "ax.set_extent ([lonW,lonE,latS,latN])\n", "ax.set_facecolor(cfeature.COLORS['water'])\n", "ax.add_feature (cfeature.STATES.with_scale(res))\n", "ax.add_feature (cfeature.RIVERS.with_scale(res))\n", "ax.add_feature (cfeature.LAND.with_scale(res))\n", "ax.add_feature (cfeature.COASTLINE.with_scale(res))\n", "ax.add_feature (cfeature.LAKES.with_scale(res))\n", "ax.add_feature (cfeature.STATES.with_scale(res))\n", "\n", "ax.add_feature(USCOUNTIES.with_scale(resCounty), linewidth=0.8, edgecolor='darkgray')\n", "\n", "# Create a station plot pointing to an Axes to draw on as well as the location of points\n", "stationplot = StationPlot(ax, lons, lats, transform=ccrs.PlateCarree(),\n", " fontsize=8)\n", "\n", "stationplot.plot_parameter('NW', tmpf, color='red')\n", "stationplot.plot_parameter('SW', rh, color='green')\n", "stationplot.plot_parameter('NE', pres, color='purple')\n", "stationplot.plot_barb(u, v,zorder=2) # zorder value set so wind barbs will display over lake features\n", "\n", "plotTitle = f'NYSM Temperature(°F), RH (%), Station Pressure (hPa), Peak 5-min Wind (kts), {titleString}'\n", "ax.set_title (plotTitle);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Tip: In the code cell above, the plotTitle string object, uses f-strings, which is a newer and better way of constructing more complex string objects in Python.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if we wanted to plot sea-level pressure (SLP) instead of station pressure? In this case, we can apply what's called a **reduction to sea-level pressure** formula. This formula requires station elevation (accounting for sensor height) in meters, temperature in Kelvin, and station pressure in hectopascals. We assume each NYSM station has its sensor height .5 meters above ground level." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Tip: At present, MetPy does not yet have a function that reduces station pressure to SLP, so we will do a units-unaware calculation. We can strip the units ... i.e., take only the magnitude of the units-aware variables via their .m attribute.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "pres.m" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "elev = nysm_data['elevation']\n", "sensorHeight = .5\n", "# Reduce station pressure to SLP. Source: https://www.sandhurstweather.org.uk/barometric.pdf \n", "slp = pres.m/np.exp(-1*(elev+sensorHeight)/((tmpc.m + 273.15) * 29.263))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "slp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a new map, substituting SLP for station pressure. We will also use the convention of the three least-significant digits to represent SLP in hectopascals." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Examples: 1018.4 hPa would be plotted as 184, while 977.2 hPa would be plotted as 772.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(18,12),dpi=150) # Increase the dots per inch from default 100 to make plot easier to read\n", "ax = fig.add_subplot(1,1,1,projection=proj)\n", "ax.set_extent ([lonW,lonE,latS,latN])\n", "ax.set_facecolor(cfeature.COLORS['water'])\n", "ax.add_feature (cfeature.STATES.with_scale(res))\n", "ax.add_feature (cfeature.RIVERS.with_scale(res))\n", "ax.add_feature (cfeature.LAND.with_scale(res))\n", "ax.add_feature (cfeature.COASTLINE.with_scale(res))\n", "ax.add_feature (cfeature.LAKES.with_scale(res))\n", "ax.add_feature (cfeature.STATES.with_scale(res))\n", "\n", "ax.add_feature(USCOUNTIES.with_scale(resCounty), linewidth=0.8, edgecolor='darkgray')\n", "\n", "stationplot = StationPlot(ax, lons, lats, transform=ccrs.PlateCarree(),\n", " fontsize=8)\n", "\n", "stationplot.plot_parameter('NW', tmpf, color='red')\n", "stationplot.plot_parameter('SW', rh, color='green')\n", "# A more complex example uses a custom formatter to control how the sea-level pressure\n", "# values are plotted. This uses the standard trailing 3-digits of the pressure value\n", "# in tenths of millibars.\n", "stationplot.plot_parameter('NE', slp, color='purple', formatter=lambda v: format(10 * v, '.0f')[-3:])\n", "stationplot.plot_barb(u, v,zorder=2)\n", "plotTitle = f'NYSM Temperature(°F), RH (%), SLP(hPa), Peak 5-min Wind (kts), {titleString}'\n", "ax.set_title (plotTitle);\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One last thing to do ... plot dewpoint instead of RH. MetPy's [dewpoint_from_relative_humidity](https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.dewpoint_from_relative_humidity.html) takes care of this!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "dwpc = dewpoint_from_relative_humidity(tmpf, rh)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "dwpc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The dewpoint is returned in units of degrees Celsius, so convert to Fahrenheit." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "dwpf = dwpc.to('degF')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the map" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = '10m'\n", "fig = plt.figure(figsize=(18,12),dpi=150) # Increase the dots per inch from default 100 to make plot easier to read\n", "ax = fig.add_subplot(1,1,1,projection=proj)\n", "ax.set_extent ([lonW,lonE,latS,latN])\n", "ax.set_facecolor(cfeature.COLORS['water'])\n", "ax.add_feature(cfeature.STATES.with_scale(res))\n", "ax.add_feature(cfeature.RIVERS.with_scale(res))\n", "ax.add_feature (cfeature.LAND.with_scale(res))\n", "ax.add_feature(cfeature.COASTLINE.with_scale(res))\n", "ax.add_feature (cfeature.LAKES.with_scale(res))\n", "ax.add_feature (cfeature.STATES.with_scale(res))\n", "\n", "ax.add_feature(USCOUNTIES.with_scale(resCounty), linewidth=0.8, edgecolor='darkgray')\n", "\n", "stationplot = StationPlot(ax, lons, lats, transform=ccrs.PlateCarree(),\n", " fontsize=8)\n", "\n", "stationplot.plot_parameter('NW', tmpf, color='red')\n", "stationplot.plot_parameter('SW', dwpf, color='green')\n", "# A more complex example uses a custom formatter to control how the sea-level pressure\n", "# values are plotted. This uses the standard trailing 3-digits of the pressure value\n", "# in tenths of millibars.\n", "stationplot.plot_parameter('NE', slp, color='purple', formatter=lambda v: format(10 * v, '.0f')[-3:])\n", "stationplot.plot_barb(u, v,zorder=2)\n", "plotTitle = f'NYSM Temperature and Dewpoint(°F), SLP (hPa), Peak 5-min Wind (kts), {titleString}'\n", "ax.set_title (plotTitle);\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Save the plot to the current directory." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "figName = f'NYSM_{figString}.png'\n", "figName" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "fig.savefig(figName)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Summary\n", "* The **MetPy** library provides methods to assign physical units to numerical arrays and perform units-aware calculations\n", "* **MetPy**'s `StationPlot` method offers a customized use of **Matplotlib**'s **Pyplot** library to plot several meteorologically-relevant parameters centered about several geo-referenced points.\n", "### What's Next?\n", "In the next notebook, we will plot METAR observations from sites across the world, using a different access method to retrieve the data.\n", "## Resources and References\n", "1. [MetPy's *calc* library](https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.html)\n", "1. [MetPy's *units* library](https://unidata.github.io/MetPy/latest/api/generated/metpy.units.html)\n", "1. [Sea-level pressure reduction formula (source: Sandhurst Weather site)](https://www.sandhurstweather.org.uk/barometric.pdf)\n", "1. [MetPy's *StationPlot* class](https://unidata.github.io/MetPy/latest/api/generated/metpy.plots.StationPlot.html#metpy.plots.StationPlot)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 Jan. 2024 Environment", "language": "python", "name": "jan24" }, "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.11.7" } }, "nbformat": 4, "nbformat_minor": 4 }