MetPy Intro: NYS Mesonet Map
Contents
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
ERROR 1: PROJ: proj_create_from_database: Open of /knight/anaconda_aug22/envs/aug22_env/share/proj failed
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.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')
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
Seriesobjects for the variables of interestExtract the underlying Numpy array via the
Series’valuesattributeAttach units to these arrays using MetPy’s
unitsclass
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']
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.
Examine these two units aware Series
wspd
| Magnitude | [5.4 4.6 3.5 8.1 6.3 2.1 5.2 6.2 4.3 7.8 5.3 8.8 3.3 9.6 3.8 0.8 3.6 4.7 |
|---|---|
| Units | meter/second |
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!
drct
| Magnitude | [359 340 304 347 320 303 358 309 326 7 307 352 353 11 308 81 297 330 |
|---|---|
| Units | degree |
Convert wind speed from m/s to knots
wspk = wspd.to('knots')
Perform the vector decomposition¶
u, v = wind_components(wspk, drct)
Take a look at one of the output components:
u
| Magnitude | [0.1831937263179225 3.058236270773797 5.640320418031061 3.541885931245742 |
|---|---|
| Units | knot |
Exercise: Create a 2-meter temperature object; attach units of degrees Celsius (degC) to it, and then convert to Fahrenheit (degF).
# Write your code here
# %load /spare11/atm350/common/mar09/mar09.py
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!¶
# 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
res = '50m'
proj = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)
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)
ax.set_extent ([lonW,lonE,latS,latN])
ax.set_facecolor(cfeature.COLORS['water'])
ax.add_feature (cfeature.STATES.with_scale(res))
ax.add_feature (cfeature.RIVERS.with_scale(res))
ax.add_feature (cfeature.LAND.with_scale(res))
ax.add_feature (cfeature.COASTLINE.with_scale(res))
ax.add_feature (cfeature.LAKES.with_scale(res))
ax.add_feature (cfeature.STATES.with_scale(res))
# Create a station plot pointing to an Axes to draw on as well as the location of points
stationplot = StationPlot(ax, lons, lats, transform=ccrs.PlateCarree(),
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($^\circ$F), RH (%), Station Pressure (hPa), Peak 5-min Wind (kts), {titleString}'
ax.set_title (plotTitle);
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.
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.
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.
pres.m
array([ 962.53, 959.89, 992.48, 1009.73, 966.7 , 1007.9 , 973.94,
976.76, 959.14, 1015.56, 999.75, 997.55, 1005.45, 1012.82,
962.38, 1007.74, 1004.69, 955.63, 1015.54, 1006.72, 1015.66,
1017.89, 981.59, 972.85, 951.1 , 1002.14, 970.62, 982.83,
952.41, 970.63, 995.47, 994.78, 973.36, 1009.03, 949.92,
985.72, 1004.91, 971.41, 979.78, 982.05, 999.7 , 979.73,
987.82, 984.07, 1018.01, 1003.89, 997.31, 958.91, 1010.65,
993.69, 977.2 , 951.31, 1011.63, 959.67, 964.49, 936.83,
997.22, 997.35, 962.74, 994.36, 1011.33, 1012.79, 963.41,
1015.23, 996.14, 1007.99, 1006.96, 975.54, 977.57, 967.1 ,
962.61, 986.28, 959.13, 968.83, 1014.09, 978.77, 983.48,
1014.88, 992.15, 969.86, 994.52, 1009.44, 957.23, 1009.78,
1012.85, 970.42, 958.26, 978.04, 1013.35, 952.77, 1003.72,
985.49, 982.42, 1010.22, 1008.89, 1018.91, 977.07, 981. ,
997.61, 1017.37, 998.84, 968.27, 1015.5 , 985.97, 1013.88,
998.36, 935.79, 1013.85, 978.87, 962.88, 980.05, 1008.13,
1007.72, 955.38, 1018.28, 959.81, 1000.32, 1008.13, 968.72,
1014.9 , 1000.42, 948.86, 970.29, 1019.04, 1011.59, 1004.2 ])
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 1025.205766
1 1024.353426
2 1027.414943
3 1020.962335
4 1025.713979
...
121 1025.108480
122 1026.126808
123 1023.728015
124 1027.146649
125 1026.637546
Name: elevation, Length: 126, 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.
Examples: 1018.4 hPa would be plotted as 184, while 977.2 hPa would be plotted as 772.
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)
ax.set_extent ([lonW,lonE,latS,latN])
ax.set_facecolor(cfeature.COLORS['water'])
ax.add_feature (cfeature.STATES.with_scale(res))
ax.add_feature (cfeature.RIVERS.with_scale(res))
ax.add_feature (cfeature.LAND.with_scale(res))
ax.add_feature (cfeature.COASTLINE.with_scale(res))
ax.add_feature (cfeature.LAKES.with_scale(res))
ax.add_feature (cfeature.STATES.with_scale(res))
stationplot = StationPlot(ax, lons, lats, transform=ccrs.PlateCarree(),
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($^\circ$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
| Magnitude | [-3.8973443699255768 -5.315778121487085 -5.300826337218609 |
|---|---|
| Units | degree_Celsius |
The dewpoint is returned in units of degrees Celsius, so convert to Fahrenheit.
dwpf = dwpc.to('degF')
Plot the map
res = '50m'
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)
ax.set_extent ([lonW,lonE,latS,latN])
ax.set_facecolor(cfeature.COLORS['water'])
ax.add_feature(cfeature.STATES.with_scale(res))
ax.add_feature(cfeature.RIVERS.with_scale(res))
ax.add_feature (cfeature.LAND.with_scale(res))
ax.add_feature(cfeature.COASTLINE.with_scale(res))
ax.add_feature (cfeature.LAKES.with_scale(res))
ax.add_feature (cfeature.STATES.with_scale(res))
stationplot = StationPlot(ax, lons, lats, transform=ccrs.PlateCarree(),
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($^\circ$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_20230309_2135.png'
fig.savefig(figName)
Summary¶
The MetPy library provides methods to assign physical units to numerical arrays and perform units-aware calculations
MetPy’s
StationPlotmethod 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, using a different access method to retrieve the data.