Introduction to MetPy

Unidata Logo

Introduction to MetPy

Unidata Python Workshop


Questions

  1. What is MetPy?
  2. How is MetPy structured?
  3. How are units handled in MetPy?

Objectives

  1. What is MetPy?
  2. Units and MetPy
  3. MetPy Constants
  4. MetPy Calculations

What is MetPy?

MetPy is a modern meteorological open-source toolkit for Python. It is a maintained project of Unidata to serve the academic meteorological community. MetPy consists of three major areas of functionality:

Plots

As meteorologists, we have many field specific plots that we make. Some of these, such as the Skew-T Log-p require non-standard axes and are difficult to plot in most plotting software. In MetPy we've baked in a lot of this specialized functionality to help you get your plots made and get back to doing science. We will go over making different kinds of plots during the workshop.

Calculations

Meteorology also has a common set of calculations that everyone ends up programming themselves. This is error-prone and a huge duplication of work! MetPy contains a set of well tested calculations that is continually growing in an effort to be at feature parity with other legacy packages such as GEMPAK.

File I/O

Finally, there are a number of odd file formats in the meteorological community. MetPy has incorporated a set of readers to help you deal with file formats that you may encounter during your research.

Units and MetPy

In order for us to discuss any of the functionality of MetPy, we first need to understand how units are inherently a part of MetPy and how to use them within this library.

Early in our scientific careers we all learn about the importance of paying attention to units in our calculations. Unit conversions can still get the best of us and have caused more than one major technical disaster, including the crash and complete loss of the $327 million Mars Climate Orbiter.

In MetPy, we use the pint library and a custom unit registry to help prevent unit mistakes in calculations. That means that every quantity you pass to MetPy should have units attached, just like if you were doing the calculation on paper! Attaching units is easy:

In [1]:
# Import the MetPy unit registry
from metpy.units import units
In [2]:
length = 10.4 * units.inches
width = 20 * units.meters
print(length, width)
10.4 inch 20 meter

Don't forget that you can use tab completion to see what units are available! Just about every imaginable quantity is there, but if you find one that isn't, we're happy to talk about adding it.

While it may seem like a lot of trouble, let's compute the area of a rectangle defined by our length and width variables above. Without units attached, you'd need to remember to perform a unit conversion before multiplying or you would end up with an area in inch-meters and likely forget about it. With units attached, the units are tracked for you.

In [3]:
area = length * width
print(area)
208.0 inch * meter

That's great, now we have an area, but it is not in a very useful unit still. Units can be converted using the .to() method. While you won't see m$^2$ in the units list, we can parse complex/compound units as strings:

In [4]:
area.to('m^2')
Out[4]:
5.2832 meter2
EXERCISE:
  • Create a variable named speed with a value of 25 knots.
  • Create a variable named time with a value of 1 fortnight.
  • Calculate how many furlongs you would travel in time at speed.
In [5]:
# YOUR CODE GOES HERE
SOLUTION
In [6]:
# %load solutions/distance.py

# Cell content replaced by load magic replacement.
speed = 25 * units.knots
time = 1 * units.fortnight
distance = speed * time
print(distance.to('furlongs'))
77332.22424242426 furlong

Temperature

Temperature units are actually relatively tricky (more like absolutely tricky as you'll see). Temperature is a non-multiplicative unit - they are in a system with a reference point. That means that not only is there a scaling factor, but also an offset. This makes the math and unit book-keeping a little more complex. Imagine adding 10 degrees Celsius to 100 degrees Celsius. Is the answer 110 degrees Celsius or 383.15 degrees Celsius (283.15 K + 373.15 K)? That's why there are delta degrees units in the unit registry for offset units. For more examples and explanation you can watch MetPy Monday #13.

Let's take a look at how this works and fails:

We would expect this to fail because we cannot add two offset units (and it does fail as an "Ambiguous operation with offset unit").

10 * units.degC + 5 * units.degC

On the other hand, we can subtract two offset quantities and get a delta:

In [7]:
10 * units.degC - 5 * units.degC
Out[7]:
5 delta_degree_Celsius

We can add a delta to an offset unit as well:

In [8]:
25 * units.degC + 5 * units.delta_degF
Out[8]:
27.77777777777778 degree_Celsius

Absolute temperature scales like Kelvin and Rankine do not have an offset and therefore can be used in addition/subtraction without the need for a delta verion of the unit.

In [9]:
273 * units.kelvin + 10 * units.kelvin
Out[9]:
283 kelvin
In [10]:
273 * units.kelvin - 10 * units.kelvin
Out[10]:
263 kelvin

Example

Let's say we're given a 12 UTC sounding, but want to know how the profile has changed when we have had several hours of diurnal heating. How do we update the surface temperature?

In [11]:
# 12 UTC temperature
temp_initial = 20 * units.degC
temp_initial
Out[11]:
20 degree_Celsius

Maybe the surface temperature increased by 5 degrees Celsius so far today - is this a temperature of 5 degC, or a temperature change of 5 degC? We subconsciously know that its a delta of 5 degC, but often write it as just adding two temperatures together, when it really is: temperature + delta(temperature)

In [12]:
# New 18 UTC temperature
temp_new = temp_initial + 5 * units.delta_degC
temp_new
Out[12]:
25 degree_Celsius
EXERCISE: A cold front is moving through, decreasing the ambient temperature of 25 degC at a rate of 2.3 degF every 10 minutes. What is the temperature after 1.5 hours?
In [13]:
# YOUR CODE GOES HERE
SOLUTION
In [14]:
# %load solutions/temperature_change.py

# Cell content replaced by load magic replacement.
temperature_change_rate = -2.3 * units.delta_degF / (10 * units.minutes)
temperature = 25 * units.degC
dt = 1.5 * units.hours
print(temperature + temperature_change_rate * dt)
13.5 degree_Celsius

Top


MetPy Constants

Another common place that problems creep into scientific code is the value of constants. Can you reproduce someone else's computations from their paper? Probably not unless you know the value of all of their constants. Was the radius of the earth 6000 km, 6300km, 6371 km, or was it actually latitude dependent?

MetPy has a set of constants that can be easily accessed and make your calculations reproducible. You can view a full table in the docs, look at the module docstring with metpy.constants? or checkout what's available with tab completion.

In [15]:
import metpy.constants as mpconst
In [16]:
mpconst.earth_avg_radius
Out[16]:
6371200.0 meter
In [17]:
mpconst.dry_air_molecular_weight
Out[17]:
28.9644 gram/mole

You may also notice in the table that most constants have a short name as well that can be used:

In [18]:
mpconst.Re
Out[18]:
6371200.0 meter
In [19]:
mpconst.Md
Out[19]:
28.9644 gram/mole

Top


MetPy Calculations

MetPy also encompasses a set of calculations that are common in meteorology (with the goal of have all of the functionality of legacy software like GEMPAK and more). The calculations documentation has a complete list of the calculations in MetPy.

We'll scratch the surface and show off a few simple calculations here, but will be using many during the workshop.

In [20]:
import metpy.calc as mpcalc
import numpy as np
In [21]:
# Make some fake data for us to work with
np.random.seed(19990503)  # So we all have the same data
u = np.random.randint(0, 15, 10) * units('m/s')
v = np.random.randint(0, 15, 10) * units('m/s')

print(u)
print(v)
[14.0 2.0 12.0 5.0 3.0 5.0 14.0 8.0 9.0 10.0] meter / second
[6.0 10.0 7.0 11.0 10.0 13.0 2.0 3.0 5.0 0.0] meter / second

Let's use the wind_direction function from MetPy to calculate wind direction from these values. Remember you can look at the docstring or the website for help.

In [22]:
direction = mpcalc.wind_direction(u, v)
print(direction)
[246.80140948635182 191.30993247402023 239.74356283647072 204.44395478041653 196.69924423399362 201.03751102542182 261.86989764584405 249.44395478041653 240.94539590092285 270.0] degree
EXERCISE:
  • Calculate the wind speed using the wind_speed function.
  • Print the wind speed in m/s and mph.
In [23]:
# YOUR CODE GOES HERE
SOLUTION
In [24]:
# %load solutions/wind_speed.py

# Cell content replaced by load magic replacement.
speed = mpcalc.wind_speed(u, v)
print(speed)
print(speed.to('mph'))
[15.231546211727817 10.198039027185569 13.892443989449804 12.083045973594572 10.44030650891055 13.92838827718412 14.142135623730951 8.54400374531753 10.295630140987 10.0] meter / second
[34.0719985051177 22.81236360769857 31.076512145333307 27.029004056895516 23.354300529953807 31.156917227058244 31.63505642387918 19.11239205734952 23.030668711943 22.36936292054402] mile_per_hour

As one final demonstration, we will calculation the dewpoint given the temperature and relative humidity:

In [25]:
mpcalc.dewpoint_from_relative_humidity(25 * units.degC, 75 * units.percent)
Out[25]:
20.264799097790046 degree_Celsius

Top