MetPy¶

Tools that you either love or hate

by Siyi Wang, 2023 Nov 20th¶

1. Introduction¶

• Meteorological Calculation: Built-in functions that calculate special physical quantities (q, lapse_rate) from basic physical quantities (T, RH, P, Wind) using well-documented equations. It allows to attach units before calculation.

• Visualize meteorology variables: Pressure, Temperature, Dew Point, Wind Direction and Wind Speed etc..

• Compatible with data in various formats: .nc, .csv, .hdf, .las from multiple data sources: Models, Observations, LIDAR, Satellite

• Built-in functions to convert units painlessly: C-F-K, sigma-mb/hpa

• Play around projections easily: Regional, Global, Polar, Half Hemisphere.

MetPy Official User Guide

MetPy Official Tutorial-MetPy Monday

2. Installing packages¶

Before installation, make sure:

MetPy currently supports Python >= 3.8, by 11/20/2023 Need the following dependencies: - "matplotlib>=3.3.0", - "numpy>=1.18.0", - "pandas>=1.0.0", - "pint>=0.15", - "pooch>=1.2.0", - "pyproj>=2.6.1", - "scipy>=1.4.0", - "traitlets>=5.0.5", - "xarray>=0.18.0"

Install using conda/mamba¶

mamba install metpyconda install metpy

Intall using conda-forge¶

conda install -c conda-forge metpy

Import MetPy and other necessary packages¶

In [124]:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

import numpy as np
from metpy.plots import SkewT
from metpy.units import units
from metpy.calc import dewpoint_from_relative_humidity ## old_name: dewpoint_rh

from metpy.calc import reduce_point_density
from metpy.plots import add_metpy_logo, current_weather, sky_cover, StationPlot

# Set up matplotlib to display plots in the notebook
%matplotlib inline

import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np from metpy.plots import SkewT from metpy.units import units from metpy.calc import dewpoint_from_relative_humidity ## old_name: dewpoint_rh from metpy.calc import reduce_point_density from metpy.plots import add_metpy_logo, current_weather, sky_cover, StationPlot # Set up matplotlib to display plots in the notebook %matplotlib inline

Using MetPy¶

Units Syntax¶

In [17]:

# First way to attach unit-multiply the units

# First way to attach unit-multiply the units

In [18]:

distance = np.arange(1, 5)
distance

distance = np.arange(1, 5) distance

Out[18]:

array([1, 2, 3, 4])

In [19]:

distance = np.arange(1, 5) * units.meters
distance

distance = np.arange(1, 5) * units.meters distance

Out[19]:

Magnitude	[1 2 3 4]
Units	meter

In [20]:

# Another way to attach units is to create the array directly with the pint.Quantity() object:

# Another way to attach units is to create the array directly with the pint.Quantity() object:

In [22]:

time = units.Quantity(np.arange(2, 10, 2), 'sec')
time

time = units.Quantity(np.arange(2, 10, 2), 'sec') time

Out[22]:

Magnitude	[2 4 6 8]
Units	second

Unit-aware Calculations¶

In [25]:

velocity = distance/time
velocity

velocity = distance/time velocity

Out[25]:

Magnitude	[0.5 0.5 0.5 0.5]
Units	meter/second

Create a Skew-T log-P plot¶

Visualizing the first type of special meteorological data, sounding data

Let's start by creating some sample data. You can replace this with your own data later.¶

In [109]:

# Pressure and Temperature data
pressure_levels = np.array([1000, 850, 700, 600, 500, 300, 200]) * units.hPa
temperature = np.array([30, 20, 10, 0, -10, -30, -40]) * units.degC
dewpoint = dewpoint_from_relative_humidity(temperature, np.array([70, 100, 80, 60, 100, 50, 40]) * units.percent)

# Pressure and Temperature data pressure_levels = np.array([1000, 850, 700, 600, 500, 300, 200]) * units.hPa temperature = np.array([30, 20, 10, 0, -10, -30, -40]) * units.degC dewpoint = dewpoint_from_relative_humidity(temperature, np.array([70, 100, 80, 60, 100, 50, 40]) * units.percent)

In [110]:

# Convert temperature and dewpoint to Kelvin for MetPy
temperature = temperature.to(units.K)
dewpoint = dewpoint.to(units.K)

# Convert temperature and dewpoint to Kelvin for MetPy temperature = temperature.to(units.K) dewpoint = dewpoint.to(units.K)

In [116]:

# Define Wind
u = np.linspace(-10, 10, len(pressure_levels)) * units.knots
v = np.array([-20, -20, -30, 5, 10, 20, 20]) * units.knots

# Define Wind u = np.linspace(-10, 10, len(pressure_levels)) * units.knots v = np.array([-20, -20, -30, 5, 10, 20, 20]) * units.knots

Plotting using SkewT() function¶

In [119]:

# Create a skew-T log-P plot
fig = plt.figure(figsize=(12, 9))
skew = SkewT(fig, rotation=45)

# Plot temperature and dewpoint profiles
skew.plot(pressure_levels, temperature, label='Temperature', color='r', linewidth=2.5)
skew.plot(pressure_levels, dewpoint, label='Dewpoint', color='g', linewidth=2.5) # linestyle='dashed'

# Add wind barbs
skew.plot_barbs(pressure_levels, u, v)

# Add other plot features
skew.ax.set_xlim(-40, 40)
skew.ax.set_ylim(1000, 100)
skew.ax.set_title('Skew-T Log-P Diagram')
skew.ax.set_xlabel('Temperature (°C)')
skew.ax.set_ylabel('Pressure (hPa)')

# Add the relevant special lines
# skew.plot_dry_adiabats()
# skew.plot_moist_adiabats()
# skew.plot_mixing_lines()

# Set x, y limits
skew.ax.set_xlim(-30, 40)
skew.ax.set_ylim(1000, 100)

# Add legend
skew.ax.legend()

# Show the plot
plt.show()

# Create a skew-T log-P plot fig = plt.figure(figsize=(12, 9)) skew = SkewT(fig, rotation=45) # Plot temperature and dewpoint profiles skew.plot(pressure_levels, temperature, label='Temperature', color='r', linewidth=2.5) skew.plot(pressure_levels, dewpoint, label='Dewpoint', color='g', linewidth=2.5) # linestyle='dashed' # Add wind barbs skew.plot_barbs(pressure_levels, u, v) # Add other plot features skew.ax.set_xlim(-40, 40) skew.ax.set_ylim(1000, 100) skew.ax.set_title('Skew-T Log-P Diagram') skew.ax.set_xlabel('Temperature (°C)') skew.ax.set_ylabel('Pressure (hPa)') # Add the relevant special lines # skew.plot_dry_adiabats() # skew.plot_moist_adiabats() # skew.plot_mixing_lines() # Set x, y limits skew.ax.set_xlim(-30, 40) skew.ax.set_ylim(1000, 100) # Add legend skew.ax.legend() # Show the plot plt.show()

For more complex graphs - adding atmospheric profile and special relevant lines¶

In [120]:

import metpy.calc as mpcalc

fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig)

# Create arrays of pressure, temperature, dewpoint, and wind components
p = [902, 897, 893, 889, 883, 874, 866, 857, 849, 841, 833, 824, 812, 796, 776, 751,
     727, 704, 680, 656, 629, 597, 565, 533, 501, 468, 435, 401, 366, 331, 295, 258,
     220, 182, 144, 106] * units.hPa
t = [-3, -3.7, -4.1, -4.5, -5.1, -5.8, -6.5, -7.2, -7.9, -8.6, -8.9, -7.6, -6, -5.1,
     -5.2, -5.6, -5.4, -4.9, -5.2, -6.3, -8.4, -11.5, -14.9, -18.4, -21.9, -25.4,
     -28, -32, -37, -43, -49, -54, -56, -57, -58, -60] * units.degC
td = [-22, -22.1, -22.2, -22.3, -22.4, -22.5, -22.6, -22.7, -22.8, -22.9, -22.4,
      -21.6, -21.6, -21.9, -23.6, -27.1, -31, -38, -44, -46, -43, -37, -34, -36,
      -42, -46, -49, -48, -47, -49, -55, -63, -72, -88, -93, -92] * units.degC
u = np.linspace(-10, 10, len(p)) * units.knots
v = np.linspace(-20, 20, len(p)) * units.knots

# Calculate parcel profile
prof = mpcalc.parcel_profile(p, t[0], td[0]).to('degC')

skew.plot(p, t, 'r')
skew.plot(p, td, 'g')
skew.plot(p, prof, 'k')  # Plot parcel profile
skew.plot_barbs(p[::5], u[::5], v[::5])

skew.ax.set_xlim(-70, 15)
skew.ax.set_ylim(1000, 100)

# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()

plt.show()

import metpy.calc as mpcalc fig = plt.figure(figsize=(9, 9)) skew = SkewT(fig) # Create arrays of pressure, temperature, dewpoint, and wind components p = [902, 897, 893, 889, 883, 874, 866, 857, 849, 841, 833, 824, 812, 796, 776, 751, 727, 704, 680, 656, 629, 597, 565, 533, 501, 468, 435, 401, 366, 331, 295, 258, 220, 182, 144, 106] * units.hPa t = [-3, -3.7, -4.1, -4.5, -5.1, -5.8, -6.5, -7.2, -7.9, -8.6, -8.9, -7.6, -6, -5.1, -5.2, -5.6, -5.4, -4.9, -5.2, -6.3, -8.4, -11.5, -14.9, -18.4, -21.9, -25.4, -28, -32, -37, -43, -49, -54, -56, -57, -58, -60] * units.degC td = [-22, -22.1, -22.2, -22.3, -22.4, -22.5, -22.6, -22.7, -22.8, -22.9, -22.4, -21.6, -21.6, -21.9, -23.6, -27.1, -31, -38, -44, -46, -43, -37, -34, -36, -42, -46, -49, -48, -47, -49, -55, -63, -72, -88, -93, -92] * units.degC u = np.linspace(-10, 10, len(p)) * units.knots v = np.linspace(-20, 20, len(p)) * units.knots # Calculate parcel profile prof = mpcalc.parcel_profile(p, t[0], td[0]).to('degC') skew.plot(p, t, 'r') skew.plot(p, td, 'g') skew.plot(p, prof, 'k') # Plot parcel profile skew.plot_barbs(p[::5], u[::5], v[::5]) skew.ax.set_xlim(-70, 15) skew.ax.set_ylim(1000, 100) # Add the relevant special lines skew.plot_dry_adiabats() skew.plot_moist_adiabats() skew.plot_mixing_lines() plt.show()

Create a Station Plot¶

Visualizing the second type of special meteorological data, on-site observation

Features: Scattered, not regularly distributed grid data

Get Data¶

In [126]:

# Import tutorial data packages
from metpy.cbook import get_test_data
from metpy.io import metar

# Import tutorial data packages from metpy.cbook import get_test_data from metpy.io import metar

In [125]:

data = metar.parse_metar_file(get_test_data('metar_20190701_1200.txt', as_file_obj=False))

# Drop rows with missing winds
data = data.dropna(how='any', subset=['wind_direction', 'wind_speed'])

data = metar.parse_metar_file(get_test_data('metar_20190701_1200.txt', as_file_obj=False)) # Drop rows with missing winds data = data.dropna(how='any', subset=['wind_direction', 'wind_speed'])

Downloading file 'metar_20190701_1200.txt' from 'https://github.com/Unidata/MetPy/raw/v1.5.1/staticdata/metar_20190701_1200.txt' to '/home/sw1225/.cache/metpy/v1.5.1'.
Downloading file 'sfstns.tbl' from 'https://github.com/Unidata/MetPy/raw/v1.5.1/staticdata/sfstns.tbl' to '/home/sw1225/.cache/metpy/v1.5.1'.
Downloading file 'master.txt' from 'https://github.com/Unidata/MetPy/raw/v1.5.1/staticdata/master.txt' to '/home/sw1225/.cache/metpy/v1.5.1'.
Downloading file 'stations.txt' from 'https://github.com/Unidata/MetPy/raw/v1.5.1/staticdata/stations.txt' to '/home/sw1225/.cache/metpy/v1.5.1'.
Downloading file 'airport-codes.csv' from 'https://github.com/Unidata/MetPy/raw/v1.5.1/staticdata/airport-codes.csv' to '/home/sw1225/.cache/metpy/v1.5.1'.

This sample data has way too many stations to plot all of them. The number of stations plotted will be reduced using reduce_point_density.

In [127]:

# Set up the map projection
proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=35,
                             standard_parallels=[35])

# Use the Cartopy map projection to transform station locations to the map and
# then refine the number of stations plotted by setting a 300km radius
point_locs = proj.transform_points(ccrs.PlateCarree(), data['longitude'].values,
                                   data['latitude'].values)
data = data[reduce_point_density(point_locs, 300000.)]

# Set up the map projection proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=35, standard_parallels=[35]) # Use the Cartopy map projection to transform station locations to the map and # then refine the number of stations plotted by setting a 300km radius point_locs = proj.transform_points(ccrs.PlateCarree(), data['longitude'].values, data['latitude'].values) data = data[reduce_point_density(point_locs, 300000.)]

Plot + reduce the number of points¶

In [129]:

# Change the DPI of the resulting figure. Higher DPI drastically improves the
# look of the text rendering.
plt.rcParams['savefig.dpi'] = 255

# Create the figure and an axes set to the projection.
fig = plt.figure(figsize=(20, 10))
add_metpy_logo(fig, 1100, 300, size='large')
ax = fig.add_subplot(1, 1, 1, projection=proj)

# Add some various map elements to the plot to make it recognizable.
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.STATES)
ax.add_feature(cfeature.BORDERS)

# Set plot bounds
ax.set_extent((-118, -73, 23, 50))

#
# Here's the actual station plot
#

# Start the station plot by specifying the axes to draw on, as well as the
# lon/lat of the stations (with transform). We also the fontsize to 12 pt.
stationplot = StationPlot(ax, data['longitude'].values, data['latitude'].values,
                          clip_on=True, transform=ccrs.PlateCarree(), fontsize=12)

# Plot the temperature and dew point to the upper and lower left, respectively, of
# the center point. Each one uses a different color.
stationplot.plot_parameter('NW', data['air_temperature'].values, color='red')
stationplot.plot_parameter('SW', data['dew_point_temperature'].values,
                           color='darkgreen')

# 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', data['air_pressure_at_sea_level'].values,
                           formatter=lambda v: format(10 * v, '.0f')[-3:])

# Plot the cloud cover symbols in the center location. This uses the codes made above and
# uses the `sky_cover` mapper to convert these values to font codes for the
# weather symbol font.
stationplot.plot_symbol('C', data['cloud_coverage'].values, sky_cover)

# Same this time, but plot current weather to the left of center, using the
# `current_weather` mapper to convert symbols to the right glyphs.
stationplot.plot_symbol('W', data['current_wx1_symbol'].values, current_weather)

# Add wind barbs
stationplot.plot_barb(data['eastward_wind'].values, data['northward_wind'].values)

# Also plot the actual text of the station id. Instead of cardinal directions,
# plot further out by specifying a location of 2 increments in x and 0 in y.
stationplot.plot_text((2, 0), data['station_id'].values)

plt.show()

# Change the DPI of the resulting figure. Higher DPI drastically improves the # look of the text rendering. plt.rcParams['savefig.dpi'] = 255 # Create the figure and an axes set to the projection. fig = plt.figure(figsize=(20, 10)) add_metpy_logo(fig, 1100, 300, size='large') ax = fig.add_subplot(1, 1, 1, projection=proj) # Add some various map elements to the plot to make it recognizable. ax.add_feature(cfeature.LAND) ax.add_feature(cfeature.OCEAN) ax.add_feature(cfeature.LAKES) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.STATES) ax.add_feature(cfeature.BORDERS) # Set plot bounds ax.set_extent((-118, -73, 23, 50)) # # Here's the actual station plot # # Start the station plot by specifying the axes to draw on, as well as the # lon/lat of the stations (with transform). We also the fontsize to 12 pt. stationplot = StationPlot(ax, data['longitude'].values, data['latitude'].values, clip_on=True, transform=ccrs.PlateCarree(), fontsize=12) # Plot the temperature and dew point to the upper and lower left, respectively, of # the center point. Each one uses a different color. stationplot.plot_parameter('NW', data['air_temperature'].values, color='red') stationplot.plot_parameter('SW', data['dew_point_temperature'].values, color='darkgreen') # 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', data['air_pressure_at_sea_level'].values, formatter=lambda v: format(10 * v, '.0f')[-3:]) # Plot the cloud cover symbols in the center location. This uses the codes made above and # uses the `sky_cover` mapper to convert these values to font codes for the # weather symbol font. stationplot.plot_symbol('C', data['cloud_coverage'].values, sky_cover) # Same this time, but plot current weather to the left of center, using the # `current_weather` mapper to convert symbols to the right glyphs. stationplot.plot_symbol('W', data['current_wx1_symbol'].values, current_weather) # Add wind barbs stationplot.plot_barb(data['eastward_wind'].values, data['northward_wind'].values) # Also plot the actual text of the station id. Instead of cardinal directions, # plot further out by specifying a location of 2 increments in x and 0 in y. stationplot.plot_text((2, 0), data['station_id'].values) plt.show()

In [ ]: