MetPy Intro: NYS Mesonet Map

Overview¶

In this notebook, we’ll use Cartopy, Matplotlib, and Pandas (with a bunch of help from MetPy) to read in, manipulate, and visualize current data from the New York State Mesonet.

Prerequisites¶

Concepts	Importance	Notes
Matplotlib	Necessary
Cartopy	Necessary
Pandas	Necessary
MetPy	Necessary	Intro

Time to learn: 30 minutes

Imports¶

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from datetime import datetime
from cartopy import crs as ccrs
from cartopy import feature as cfeature
from metpy.calc import wind_components, dewpoint_from_relative_humidity
from metpy.units import units
from metpy.plots import StationPlot, USCOUNTIES

Create a Pandas DataFrame object pointing to the most recent set of NYSM obs.

nysm_data = pd.read_csv('https://www.atmos.albany.edu/products/nysm/nysm_latest.csv')

nysm_data

nysm_data.columns

Index(['station', 'time', 'temp_2m [degC]', 'temp_9m [degC]',
       'relative_humidity [percent]', 'precip_incremental [mm]',
       'precip_local [mm]', 'precip_max_intensity [mm/min]',
       'avg_wind_speed_prop [m/s]', 'max_wind_speed_prop [m/s]',
       'wind_speed_stddev_prop [m/s]', 'wind_direction_prop [degrees]',
       'wind_direction_stddev_prop [degrees]', 'avg_wind_speed_sonic [m/s]',
       'max_wind_speed_sonic [m/s]', 'wind_speed_stddev_sonic [m/s]',
       'wind_direction_sonic [degrees]',
       'wind_direction_stddev_sonic [degrees]', 'solar_insolation [W/m^2]',
       'station_pressure [mbar]', 'snow_depth [cm]', 'frozen_soil_05cm [bit]',
       'frozen_soil_25cm [bit]', 'frozen_soil_50cm [bit]',
       'soil_temp_05cm [degC]', 'soil_temp_25cm [degC]',
       'soil_temp_50cm [degC]', 'soil_moisture_05cm [m^3/m^3]',
       'soil_moisture_25cm [m^3/m^3]', 'soil_moisture_50cm [m^3/m^3]', 'lat',
       'lon', 'elevation', 'name'],
      dtype='object')

for col in nysm_data.columns:
    print(col)

station
time
temp_2m [degC]
temp_9m [degC]
relative_humidity [percent]
precip_incremental [mm]
precip_local [mm]
precip_max_intensity [mm/min]
avg_wind_speed_prop [m/s]
max_wind_speed_prop [m/s]
wind_speed_stddev_prop [m/s]
wind_direction_prop [degrees]
wind_direction_stddev_prop [degrees]
avg_wind_speed_sonic [m/s]
max_wind_speed_sonic [m/s]
wind_speed_stddev_sonic [m/s]
wind_direction_sonic [degrees]
wind_direction_stddev_sonic [degrees]
solar_insolation [W/m^2]
station_pressure [mbar]
snow_depth [cm]
frozen_soil_05cm [bit]
frozen_soil_25cm [bit]
frozen_soil_50cm [bit]
soil_temp_05cm [degC]
soil_temp_25cm [degC]
soil_temp_50cm [degC]
soil_moisture_05cm [m^3/m^3]
soil_moisture_25cm [m^3/m^3]
soil_moisture_50cm [m^3/m^3]
lat
lon
elevation
name

Create several Series objects for some of the columns.

stid = nysm_data['station']
lats = nysm_data['lat']
lons = nysm_data['lon']

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 method.

This method requires us to do the following:

Create Pandas Series objects for the variables of interest
Extract the underlying Numpy array via the Series’ values attribute
Attach units to these arrays using MetPy’s units class

Perform these three steps¶

wspd = nysm_data['max_wind_speed_prop [m/s]'].values * units['m/s']
drct = nysm_data['wind_direction_prop [degrees]'].values * units['degree']

Examine these two units aware Series

wspd

drct

Convert wind speed from m/s to knots

wspk = wspd.to('knots')

wspk

Perform the vector decomposition¶

u, v = wind_components(wspk, drct)

Take a look at one of the output components:

# Write your code here
tmpc = nysm_data['temp_2m [degC]'].values * units('degC')
tmpf = tmpc.to('degF')

Click to reveal only AFTER you have tried your own code!

tmpc = nysm_data['temp_2m [degC]'].values * units('degC')
tmpf = tmpc.to('degF')

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.

Create units-aware objects for relative humidity and station pressure

rh = nysm_data['relative_humidity [percent]'].values * units('percent')
pres = nysm_data['station_pressure [mbar]'].values * units('mbar')

Plot the map, centered over NYS, with add some geographic features, and the mesonet data.¶

Determine the current time, for use in the plot’s title and a version to be saved to disk.

timeString = nysm_data['time'][0]
timeObj = datetime.strptime(timeString,"%Y-%m-%d %H:%M:%S")
titleString = datetime.strftime(timeObj,"%B %d %Y, %H%M UTC")
figString = datetime.strftime(timeObj,"%Y%m%d_%H%M")

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!¶

Be patient: this may take a minute or so to plot, if you chose the highest resolution for the shapefiles!¶

Specify two resolutions: one for the Natural Earth shapefiles, and the other for MetPy’s US County shapefiles.¶

res_nearth = '10m'
res_county = '5m'

# Set the domain for defining the plot region.
latN = 45.2
latS = 40.2
lonW = -80.0
lonE = -72.0
cLat = (latN + latS)/2
cLon = (lonW + lonE )/2

proj_map = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)

# Specify the dataset's map projection
proj_data = ccrs.PlateCarree()

fig = plt.figure(figsize=(18,12),dpi=150) # Increase the dots per inch from default 100 to make plot easier to read
ax = fig.add_subplot(1,1,1,projection=proj_map)
ax.set_extent ([lonW,lonE,latS,latN])
ax.set_facecolor(cfeature.COLORS['water'])
ax.add_feature (cfeature.STATES.with_scale(res_nearth))
ax.add_feature (cfeature.RIVERS.with_scale(res_nearth))
ax.add_feature (cfeature.LAND.with_scale(res_nearth))
ax.add_feature (cfeature.COASTLINE.with_scale(res_nearth))
ax.add_feature (cfeature.LAKES.with_scale(res_nearth))

ax.add_feature(USCOUNTIES.with_scale(res_county), linewidth=0.8, edgecolor='darkgray')

# Create a station plot pointing to an Axes to draw on as well as the location of points
stationplot = StationPlot(ax, lons, lats, transform=proj_data,
                          fontsize=8)

stationplot.plot_parameter('NW', tmpf, color='red')
stationplot.plot_parameter('SW', rh, color='green')
stationplot.plot_parameter('NE', pres, color='purple')
stationplot.plot_barb(u, v,zorder=2) # zorder value set so wind barbs will display over lake features

plotTitle = f'NYSM Temperature(°F), RH (%), Station Pressure (hPa), Peak 5-min Wind (kts), {titleString}'
ax.set_title (plotTitle);

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.

pres.m

array([ 950.68,  946.9 ,  976.9 ,  995.08,  952.12,  991.69,  961.81,
        963.47,  946.55, 1002.94,  983.8 ,  985.01,  988.84,  999.73,
        949.82,  989.81,  988.35,  944.37,  997.3 ,  991.3 ,  997.92,
        999.22,  965.69,  960.11,  940.04,  986.64,  958.29,  966.84,
        941.31,  957.47,  980.95,  979.02,  959.76,  992.83,  938.67,
        969.49,  990.8 ,  957.08,  965.51,  966.86,  983.06,  966.34,
        971.47,  970.32,  997.93,  988.6 ,  981.81,  945.74,  992.03,
        979.05,  964.19,  940.05,  993.73,  946.18,  951.33,  926.85,
        981.23,  982.83,  949.49,  977.7 ,  994.91,  995.62,  948.97,
        936.9 ,  996.77,  978.96,  995.47,  990.05,  961.49,  964.47,
        954.05,  949.58,  969.92,  946.31,  957.41,  997.01,  963.65,
        970.33,  998.03,  978.52,  957.35,  980.05,  992.  ,  944.49,
        991.54, 1000.32,  958.46,  945.61,  964.68,  996.98,  940.56,
        987.83,  969.54,  968.83,  992.4 ,  992.03,  999.98,  963.81,
        967.59,  985.07, 1006.66,  980.75,  954.72, 1002.44,  971.56,
       1002.03,  985.78,  926.31,  993.55,  965.56,  949.61,  966.85,
        990.55,  992.67,  942.35, 1006.4 ,  947.21,  987.11,  992.1 ,
        956.58,  996.9 ,  985.16,  936.25,  957.04,  998.99,  994.78,
        988.32])

elev = nysm_data['elevation']
sensorHeight = .5
# Reduce station pressure to SLP. Source: https://www.sandhurstweather.org.uk/barometric.pdf 
slp = pres.m/np.exp(-1*(elev+sensorHeight)/((tmpc.m + 273.15) * 29.263))

slp

0      1009.392569
1      1005.233540
2      1009.364593
3      1005.411178
4      1005.081443
          ...     
122    1007.124861
123    1008.683298
124    1003.349826
125    1009.248986
126    1009.325232
Name: elevation, Length: 127, dtype: float64

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.

fig = plt.figure(figsize=(18,12),dpi=150) # Increase the dots per inch from default 100 to make plot easier to read
ax = fig.add_subplot(1,1,1,projection=proj_map)
ax.set_extent ([lonW,lonE,latS,latN])
ax.set_facecolor(cfeature.COLORS['water'])
ax.add_feature (cfeature.STATES.with_scale(res_nearth))
ax.add_feature (cfeature.RIVERS.with_scale(res_nearth))
ax.add_feature (cfeature.LAND.with_scale(res_nearth))
ax.add_feature (cfeature.COASTLINE.with_scale(res_nearth))
ax.add_feature (cfeature.LAKES.with_scale(res_nearth))

ax.add_feature(USCOUNTIES.with_scale(res_county), linewidth=0.8, edgecolor='darkgray')

stationplot = StationPlot(ax, lons, lats, transform=proj_data,
                          fontsize=8)

stationplot.plot_parameter('NW', tmpf, color='red')
stationplot.plot_parameter('SW', rh, color='green')

# A more complex example uses a custom formatter to control how the sea-level pressure
# values are plotted. This uses the standard trailing 3-digits of the pressure value
# in tenths of millibars.
stationplot.plot_parameter('NE', slp, color='purple', formatter=lambda v: format(10 * v, '.0f')[-3:])

stationplot.plot_barb(u, v,zorder=2)
plotTitle = f'NYSM Temperature(°F), RH (%), SLP(hPa), Peak 5-min Wind (kts), {titleString}'
ax.set_title (plotTitle);

One last thing to do ... plot dewpoint instead of RH. MetPy’s dewpoint_from_relative_humidity takes care of this!

dwpc = dewpoint_from_relative_humidity(tmpf, rh)

dwpc

The dewpoint is returned in units of degrees Celsius, so convert to Fahrenheit.

dwpf = dwpc.to('degF')

Plot the map

fig = plt.figure(figsize=(18,12),dpi=150) # Increase the dots per inch from default 100 to make plot easier to read
ax = fig.add_subplot(1,1,1,projection=proj_map)
ax.set_extent ([lonW,lonE,latS,latN])
ax.set_facecolor(cfeature.COLORS['water'])
ax.add_feature(cfeature.STATES.with_scale(res_nearth))
ax.add_feature(cfeature.RIVERS.with_scale(res_nearth))
ax.add_feature (cfeature.LAND.with_scale(res_nearth))
ax.add_feature(cfeature.COASTLINE.with_scale(res_nearth))
ax.add_feature (cfeature.LAKES.with_scale(res_nearth))

ax.add_feature(USCOUNTIES.with_scale(res_county), linewidth=0.8, edgecolor='darkgray')

stationplot = StationPlot(ax, lons, lats, transform=proj_data,
                          fontsize=8)

stationplot.plot_parameter('NW', tmpf, color='red')
stationplot.plot_parameter('SW', dwpf, color='green')
# A more complex example uses a custom formatter to control how the sea-level pressure
# values are plotted. This uses the standard trailing 3-digits of the pressure value
# in tenths of millibars.
stationplot.plot_parameter('NE', slp, color='purple', formatter=lambda v: format(10 * v, '.0f')[-3:])
stationplot.plot_barb(u, v,zorder=2)
plotTitle = f'NYSM Temperature and Dewpoint(°F), SLP (hPa), Peak 5-min Wind (kts), {titleString}'
ax.set_title (plotTitle);

Save the plot to the current directory.

figName = f'NYSM_{figString}.png'
figName

'NYSM_20260417_0015.png'

fig.savefig(figName)

Summary¶

The MetPy library provides methods to assign physical units to numerical arrays and perform units-aware calculations
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.

What’s Next?¶

In the next notebook, we will plot METAR observations from sites across the world.

Overview¶

Prerequisites¶

Imports¶

Perform these three steps¶

Perform the vector decomposition¶

Plot the map, centered over NYS, with add some geographic features, and the mesonet data.¶

Be patient: this may take a minute or so to plot, if you chose the highest resolution for the shapefiles!¶

Specify two resolutions: one for the Natural Earth shapefiles, and the other for MetPy’s US County shapefiles.¶

Summary¶

What’s Next?¶

Resources and References¶