Getting started with MetPy#

Before we get started, we test the learning environment and the most important packages needed to run the notebooks. This is not so much a continuous coherent exercise as individual examples based on the different packages.This exercise is neither an introduction to Python nor extensive tutorials for the individual packages. I advise you, if you have little or no experience with the packages, to work through the relevant tutorial on the websites. All packages offer very good and extensive tutorials. Most of the functions presented here have been taken from these websites.

Learning objectives:
  • Getting to know the learning environment
  • Testing the MetPy package
  • Very brief overview of the function of the packages
How to proceed:
  • Testing MetPy

MetPy is a collection of tools in Python for reading, visualizing, and performing calculations with weather data. One of the most significant differences in syntax for MetPy, compared to other Python libraries, is the frequent requirement of units to be attached to arrays before being passed to MetPy functions. There are very few exceptions to this, and you’ll usually be safer to always use units whenever applicable to make sure that your analyses are done correctly. Once you get used to the units syntax, it becomes very handy, as you never have to worry about unit conversion for any calculation. MetPy does it for you!

Let’s load the MetPy and numpy library

# Load the pandas package
import numpy as np
from metpy.units import units

Please note, we have loaded only the units module here. We can assign a unit to an array by

# Attach the unit meters to the distance array
distance = np.arange(1, 5) * units.meters

Similarly, we can attach create a time array

# This is another way to attach units to an array 
time = units.Quantity(np.arange(2, 10, 2), 'sec')

Now, we can simply do unit-aware calculations

# Calculate the velocity
distance/time
\[\begin{pmatrix}0.5 & 0.5 & 0.5 & 0.5\end{pmatrix}\ \frac{\mathrm{meter}}{\mathrm{second}}\]

MetPy with xarray#

MetPy works great with xarray. MetPy’s suite of meteorological calculations are designed to integrate with xarray DataArrays and provides DataArray and Dataset accessors (collections of methods and properties attached to the .metpy property) for coordinate/CRS and unit operations.

First, some imports …

# Import xarray
import xarray as xr

# Any import of metpy activates the accessors
import metpy
from metpy.cbook import get_test_data

We can open a netCDF-file using xarray

# Open the netCDF file as a xarray Dataset
ds = xr.open_dataset(get_test_data('irma_gfs_example.nc', as_file_obj = False))

# View a summary of the Dataset
ds

<xarray.Dataset>
Dimensions:                              (time1: 9, latitude: 81,
                                          isobaric3: 31, isobaric1: 21,
                                          longitude: 131)
Coordinates:
  * time1                                (time1) datetime64[ns] 2017-09-05T12...
    reftime                              datetime64[ns] ...
  * latitude                             (latitude) float32 50.0 49.5 ... 10.0
  * isobaric3                            (isobaric3) float64 100.0 ... 1e+05
  * isobaric1                            (isobaric1) float64 1e+04 ... 1e+05
  * longitude                            (longitude) float32 250.0 ... 315.0
Data variables:
    Vertical_velocity_pressure_isobaric  (time1, isobaric1, latitude, longitude) float32 ...
    Relative_humidity_isobaric           (time1, isobaric3, latitude, longitude) float32 ...
    Temperature_isobaric                 (time1, isobaric3, latitude, longitude) float32 ...
    u-component_of_wind_isobaric         (time1, isobaric3, latitude, longitude) float32 ...
    v-component_of_wind_isobaric         (time1, isobaric3, latitude, longitude) float32 ...
    Geopotential_height_isobaric         (time1, isobaric3, latitude, longitude) float32 ...
    LatLon_361X720-0p25S-180p00E         int32 ...
Attributes: (12/13)
    Originating_or_generating_Center:                                        ...
    Originating_or_generating_Subcenter:                                     ...
    GRIB_table_version:                                                      ...
    Type_of_generating_process:                                              ...
    Analysis_or_forecast_generating_process_identifier_defined_by_originating...
    Conventions:                                                             ...
    ...                                                                                ...
    featureType:                                                             ...
    History:                                                                 ...
    geospatial_lat_min:                                                      ...
    geospatial_lat_max:                                                      ...
    geospatial_lon_min:                                                      ...
    geospatial_lon_max:                                                      ...

This Dataset consists of dimensions and their associated coordinates, which in turn make up the axes along which the data variables are defined. The dataset also has a dictionary-like collection of attributes. What happens if we look at just a single data variable?

temperature = ds['Temperature_isobaric']
temperature

<xarray.DataArray 'Temperature_isobaric' (time1: 9, isobaric3: 31,
                                          latitude: 81, longitude: 131)>
[2960469 values with dtype=float32]
Coordinates:
  * time1      (time1) datetime64[ns] 2017-09-05T12:00:00 ... 2017-09-06T12:0...
    reftime    datetime64[ns] ...
  * latitude   (latitude) float32 50.0 49.5 49.0 48.5 ... 11.5 11.0 10.5 10.0
  * isobaric3  (isobaric3) float64 100.0 200.0 300.0 ... 9.5e+04 9.75e+04 1e+05
  * longitude  (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0
Attributes:
    long_name:                      Temperature @ Isobaric surface
    units:                          K
    Grib_Variable_Id:               VAR_0-0-0_L100
    Grib2_Parameter:                [0 0 0]
    Grib2_Parameter_Discipline:     Meteorological products
    Grib2_Parameter_Category:       Temperature
    Grib2_Parameter_Name:           Temperature
    Grib2_Level_Type:               100
    Grib2_Level_Desc:               Isobaric surface
    Grib2_Generating_Process_Type:  Forecast
    grid_mapping:                   LatLon_361X720-0p25S-180p00E

This is a DataArray, which stores just a single data variable with its associated coordinates and attributes. These individual DataArrays are the kinds of objects that MetPy’s calculations take as input (more on that in Calculations section below).

Coordinates and Corrdinate Reference Systems#

MetPy’s first set of helpers comes with identifying coordinate types. In a given dataset, coordinates can have a variety of different names (time, isobaric1 …). Following CF conventions, as well as using some fall-back regular expressions, MetPy can systematically identify coordinates of the following types:

  • time
  • vertical
  • latitude
  • y
  • longitude
  • x

When identifying a single coordinate, it is best to use the property directly associated with that type

temperature.metpy.time

<xarray.DataArray 'time1' (time1: 9)>
array(['2017-09-05T12:00:00.000000000', '2017-09-05T15:00:00.000000000',
       '2017-09-05T18:00:00.000000000', '2017-09-05T21:00:00.000000000',
       '2017-09-06T00:00:00.000000000', '2017-09-06T03:00:00.000000000',
       '2017-09-06T06:00:00.000000000', '2017-09-06T09:00:00.000000000',
       '2017-09-06T12:00:00.000000000'], dtype='datetime64[ns]')
Coordinates:
  * time1    (time1) datetime64[ns] 2017-09-05T12:00:00 ... 2017-09-06T12:00:00
    reftime  datetime64[ns] ...
Attributes:
    standard_name:  time
    long_name:      time
    udunits:        Hour since 2017-09-05T12:00:00Z
    _metpy_axis:    time
x, y = temperature.metpy.coordinates('x', 'y')
x

<xarray.DataArray 'longitude' (longitude: 131)>
array([250. , 250.5, 251. , 251.5, 252. , 252.5, 253. , 253.5, 254. , 254.5,
       255. , 255.5, 256. , 256.5, 257. , 257.5, 258. , 258.5, 259. , 259.5,
       260. , 260.5, 261. , 261.5, 262. , 262.5, 263. , 263.5, 264. , 264.5,
       265. , 265.5, 266. , 266.5, 267. , 267.5, 268. , 268.5, 269. , 269.5,
       270. , 270.5, 271. , 271.5, 272. , 272.5, 273. , 273.5, 274. , 274.5,
       275. , 275.5, 276. , 276.5, 277. , 277.5, 278. , 278.5, 279. , 279.5,
       280. , 280.5, 281. , 281.5, 282. , 282.5, 283. , 283.5, 284. , 284.5,
       285. , 285.5, 286. , 286.5, 287. , 287.5, 288. , 288.5, 289. , 289.5,
       290. , 290.5, 291. , 291.5, 292. , 292.5, 293. , 293.5, 294. , 294.5,
       295. , 295.5, 296. , 296.5, 297. , 297.5, 298. , 298.5, 299. , 299.5,
       300. , 300.5, 301. , 301.5, 302. , 302.5, 303. , 303.5, 304. , 304.5,
       305. , 305.5, 306. , 306.5, 307. , 307.5, 308. , 308.5, 309. , 309.5,
       310. , 310.5, 311. , 311.5, 312. , 312.5, 313. , 313.5, 314. , 314.5,
       315. ], dtype=float32)
Coordinates:
    reftime    datetime64[ns] ...
  * longitude  (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0
Attributes:
    units:          degrees_east
    standard_name:  longitude
    _metpy_axis:    x,longitude

These coordinate type aliases can also be used in MetPy’s wrapped .sel and .loc for indexing and selecting on DataArrays. For example, to access 500 hPa heights at 1800Z,

heights = ds['Geopotential_height_isobaric'].metpy.sel(
    time='2017-09-05 18:00',
    vertical=50000.
)
heights

<xarray.DataArray 'Geopotential_height_isobaric' (latitude: 81, longitude: 131)>
[10611 values with dtype=float32]
Coordinates:
    time1      datetime64[ns] 2017-09-05T18:00:00
    reftime    datetime64[ns] ...
  * latitude   (latitude) float32 50.0 49.5 49.0 48.5 ... 11.5 11.0 10.5 10.0
    isobaric3  float64 5e+04
  * longitude  (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0
Attributes:
    long_name:                      Geopotential height @ Isobaric surface
    units:                          gpm
    Grib_Variable_Id:               VAR_0-3-5_L100
    Grib2_Parameter:                [0 3 5]
    Grib2_Parameter_Discipline:     Meteorological products
    Grib2_Parameter_Category:       Mass
    Grib2_Parameter_Name:           Geopotential height
    Grib2_Level_Type:               100
    Grib2_Level_Desc:               Isobaric surface
    Grib2_Generating_Process_Type:  Forecast
    grid_mapping:                   LatLon_361X720-0p25S-180p00E

Beyond just the coordinates themselves, a common need for both calculations with and plots of geospatial data is knowing the coordinate reference system (CRS) on which the horizontal spatial coordinates are defined. MetPy follows the CF Conventions for its CRS definitions, which it then caches on the metpy_crs coordinate in order for it to persist through calculations and other array operations. There are two ways to do so in MetPy:

First, if your dataset is already conforming to the CF Conventions, it will have a grid mapping variable that is associated with the other data variables by the grid_mapping attribute. This is automatically parsed via the .parse_cf() method:

# Parse full dataset
data_parsed = ds.metpy.parse_cf()

# Parse subset of dataset
data_subset = ds.metpy.parse_cf([
    'u-component_of_wind_isobaric',
    'v-component_of_wind_isobaric',
    'Vertical_velocity_pressure_isobaric'
])

# Parse single variable
relative_humidity = ds.metpy.parse_cf('Relative_humidity_isobaric')

Notice the newly added metpy_crs non-dimension coordinate. Now how can we use this in practice? For individual DataArrays, we can access the cartopy and pyproj objects corresponding to this CRS:

# Cartopy CRS, useful for plotting
relative_humidity.metpy.cartopy_crs
2023-04-19T08:34:15.514125 image/svg+xml Matplotlib v3.7.0, https://matplotlib.org/ 
<cartopy.crs.PlateCarree object at 0x16dcd7f10>
# pyproj CRS, useful for projection transformations and forward/backward azimuth and great
# circle calculations
relative_humidity.metpy.pyproj_crs

<Geographic 2D CRS: {"$schema": "https://proj.org/schemas/v0.2/projjso ...>
Name: undefined
Axis Info [ellipsoidal]:
- lon[east]: Longitude (degree)
- lat[north]: Latitude (degree)
Area of Use:
- undefined
Datum: undefined
- Ellipsoid: undefined
- Prime Meridian: Greenwich

Units#

Since unit-aware calculations are a major part of the MetPy library, unit support is a major part of MetPy’s xarray integration!

One very important point of consideration is that xarray data variables (in both Datasets and DataArrays) can store both unit-aware and unit-naive array types. Unit-naive array types will be used by default in xarray, so we need to convert to a unit-aware type if we want to use xarray operations while preserving unit correctness. MetPy provides the .quantify() method for this (named since we are turning the data stored inside the xarray object into a Pint Quantity object)

temperature = data_parsed['Temperature_isobaric'].metpy.quantify()
temperature

<xarray.DataArray 'Temperature_isobaric' (time1: 9, isobaric3: 31,
                                          latitude: 81, longitude: 131)>
<Quantity([[[[258.      258.      257.9     ... 258.2     258.2     258.3    ]
   [257.9     257.9     257.8     ... 258.2     258.2     258.2    ]
   [257.8     257.8     257.8     ... 258.1     258.1     258.1    ]
   ...
   [260.8     260.8     260.8     ... 258.      258.1     258.2    ]
   [261.      261.      261.      ... 258.1     258.1     258.2    ]
   [261.2     261.2     261.2     ... 258.1     258.2     258.3    ]]

  [[256.5     256.5     256.5     ... 255.6     255.6     255.5    ]
   [256.4     256.4     256.3     ... 255.5     255.5     255.5    ]
   [256.2     256.2     256.2     ... 255.4     255.4     255.4    ]
   ...
   [259.1     259.1     259.1     ... 261.1     261.4     261.6    ]
   [259.2     259.2     259.2     ... 261.1     261.4     261.6    ]
   [259.4     259.3     259.3     ... 261.1     261.4     261.6    ]]

  [[251.4     251.4     251.4     ... 250.5     250.5     250.5    ]
   [251.2     251.3     251.3     ... 250.4     250.4     250.4    ]
   [251.1     251.1     251.1     ... 250.3     250.3     250.2    ]
   ...
...
   ...
   [295.40698 295.40698 295.40698 ... 295.70697 295.70697 295.507  ]
   [295.30698 295.30698 295.30698 ... 295.20697 295.30698 295.40698]
   [295.30698 295.107   295.107   ... 294.80698 295.007   295.107  ]]

  [[291.9     290.2     290.5     ... 288.4     288.6     288.9    ]
   [291.4     290.7     290.8     ... 289.3     289.5     289.9    ]
   [290.7     290.4     290.5     ... 289.9     290.1     290.5    ]
   ...
   [297.5     297.6     297.5     ... 297.9     297.8     297.7    ]
   [297.4     297.4     297.3     ... 297.4     297.5     297.6    ]
   [297.2     297.1     297.      ... 297.      297.2     297.3    ]]

  [[293.3     291.6     291.9     ... 286.8     286.9     287.     ]
   [292.8     292.1     292.2     ... 287.      286.9     287.     ]
   [292.1     291.8     291.9     ... 287.1     287.3     287.5    ]
   ...
   [299.6     299.7     299.7     ... 300.1     300.      299.8    ]
   [299.5     299.5     299.4     ... 299.5     299.7     299.7    ]
   [299.4     299.3     299.2     ... 299.2     299.3     299.4    ]]]], 'kelvin')>
Coordinates:
  * time1      (time1) datetime64[ns] 2017-09-05T12:00:00 ... 2017-09-06T12:0...
    reftime    datetime64[ns] 2017-09-05T12:00:00
  * latitude   (latitude) float32 50.0 49.5 49.0 48.5 ... 11.5 11.0 10.5 10.0
  * longitude  (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0
    metpy_crs  object Projection: latitude_longitude
  * isobaric3  (isobaric3) float64 100.0 200.0 300.0 ... 9.5e+04 9.75e+04 1e+05
Attributes:
    long_name:                      Temperature @ Isobaric surface
    Grib_Variable_Id:               VAR_0-0-0_L100
    Grib2_Parameter:                [0 0 0]
    Grib2_Parameter_Discipline:     Meteorological products
    Grib2_Parameter_Category:       Temperature
    Grib2_Parameter_Name:           Temperature
    Grib2_Level_Type:               100
    Grib2_Level_Desc:               Isobaric surface
    Grib2_Generating_Process_Type:  Forecast
    grid_mapping:                   LatLon_361X720-0p25S-180p00E

Notice how the units are now represented in the data itself, rather than as a text attribute. Now, even if we perform some kind of xarray operation (such as taking the zonal mean), the units are preserved.

# Take the mean over time at the 1000 hPa level
temperature.sel(isobaric3=1000).mean('time1')

Plotting#

# Load data
ds = xr.open_dataset(get_test_data('narr_example.nc', as_file_obj = False))

# Parse full dataset
ds = ds.metpy.parse_cf()

# Grab lat/lon values from file as unit arrays
lats = ds.lat.metpy.unit_array
lons = ds.lon.metpy.unit_array

# Get the valid time
vtime = ds.Temperature.metpy.time[0]

# Get the 700-hPa heights without manually identifying the vertical coordinate
hght_700 = ds.Geopotential_height.metpy.sel(vertical=700 * units.hPa,
                                                 time=vtime)

# Cartopy CRS, useful for plotting
hght_700.metpy.cartopy_crs

# Import cartopy for plotting
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

# Open figure
fig = plt.figure(figsize=(12, 12))

# Add a plot with coordinate axes (crs)
ax = fig.add_subplot(1, 1, 1, projection=hght_700.metpy.cartopy_crs)

# Get coordinates
x = ds.x
y = ds.y

# Plot data
ax.imshow(hght_700, extent=(x.min(), x.max(), y.min(), y.max()),
          cmap='RdBu', origin='lower' if y[0] < y[-1] else 'upper')

# Add coastlines
ax.coastlines(color='tab:blue', resolution='10m')

# Show the plot
plt.show()

Calculations#

import numpy as np
from metpy.units import units
import metpy.calc as mpcalc

temperature = [20] * units.degC
rel_humidity  = [50] * units.percent
print(mpcalc.dewpoint_from_relative_humidity(temperature, rel_humidity))

speed = np.array([5, 10, 15, 20]) * units.knots
direction = np.array([0, 90, 180, 270]) * units.degrees
u, v = mpcalc.wind_components(speed, direction)
print(u, v)
Reminder
  • Import the package, aka import metpy
  • Metpy works with units
  • MetPy aims to have three primary purposes: read and write meteorological data (I/O), calculate meteorological quantities with well-documented equations, and create publication-quality plots of meteorological data
Homework: Check out the metpy tutorial and get familiar with the syntax.