Exploring MetPy¶

As atmospheric scientists, you view, analyze, and synthesize a multitude of datasets every day. These data and model output are collected by numerous entities and are stored in disparate locatons. Luckily, we have access to the THREDDS catalog which nicely aggregates and organizes that data for us. Using MetPy, we can programmatically access, analyze, and display those data for our specific needs on demand.

Table of Contents ¶

• Objective
• Strategy
• Step 0: Import required packages
• Step 1: Browse the THREDDS Data Server (TDS)
• Step 2: Satellite data
  • Step 2a: Obtain Satellite data
  • Step 2b: Prepare Satellite data
  • Step 2c: Visualize Satellite data
• Step 3: Model output, equivalent potential temperature ($\theta_e$)
  • Step 3a: Obtain Model Output
  • Step 3b: Prepare Model Output
  • Step 3c: Visualize Model Output
• Step 4: Surface observations
  • Step 4a: Obtain Surface observations
  • Step 4b: Prepare Surface observations
  • Step 4c: Visualize Surface observations
• Step 5: Create multi-layer plot

Objective¶

Today, you will explore a few of the many uses of MetPy by creating a multi-layer plot like this one. Plots like these are ubiquitous in research, forecasting, and education. You may have seen a plot like this in a textbook!

Strategy ¶

Creating plots like this using MetPy follows this workflow:

Since the plot you will be creating today includes three layers (satellite imagery, model output, and surface observations), you will repeat this process three times. With each iteration, you will discover nuances of each data product, what the product looks like on its own, and how to prepare it for a final production-quality plot.

Step 0: Import required packages ¶

Top

But first, we need to import all our required packages. Today we're working with:

• datetime
• io
• cartopy
• matplotlib
• metpy
• siphon

Step 1: Browse the THREDDS Data Server (TDS) ¶

Top

Breakout room activity

The THREDDS Data Server provides us with coherent access to a large collection of real-time and archived datasets from a variety of environmental data sources at a number of distributed server sites. You can browse the TDS in your web browser using this link:

https://thredds.ucar.edu/

Instructions In your breakout rooms, take a few moments to browse the catalog in a new tab of your web browser. Then, as a group, find the following information in the cell below. Determine one representative from each room to log the information in the cell below and share with the larger group at the end of the activity.

1. Find the URL for current GOES East Cloud and Moisture Imagery, Channel 02, CONUS extent

Answer: https://thredds.ucar.edu/thredds/catalog/satellite/goes/east/products/CloudAndMoistureImagery/CONUS/Channel02/current/catalog.html

2. What data format are these data stored in?

Answer: netCDF-4

3. Find the URL for the full collection of the Real Time Mesoscale Analysis output

Answer: https://thredds.ucar.edu/thredds/catalog/grib/NCEP/RTMA/CONUS_2p5km/catalog.html?dataset=grib/NCEP/RTMA/CONUS_2p5km/TP

4. What data format are these data stored in?

Answer: GRIB-2

Step 2: Satellite Data ¶

Top

The first layer we will obtain is satellite imagery. To do this, we will create a THREDDS Data Server (TDS) catalog object, then extract a single image with its metadata, and visualize the result in a simple plot. For these tasks, we will use xarray to store the information in memory.

Step 2a: Obtain Satellite Data ¶

Top

Step 2b: Prepare Satellite Data ¶

Top

Step 2c: Visualize Satellite Data ¶

Top

ACTIVITY: Plot CMI¶

Now that we have the data ready, it's time to plot. Incrementally plotting data is considered a best practice for understanding the data you are working with and determining if you are on your way to creating your intended product. To plot this image, we will use matplotlib to display our data. Let's first make sure we have the information we need:

• What is the name of the dataarray variable we need to plot?
• What is the matplotlib function to use?
• What is the correct syntax for using the function?

Step 3: Model output, equivalent potential temperature ($\theta_e$) ¶

Top

Now we grab data from the Real-time Meso-analysis (RTMA), which gives us a gridded estimate of the realtime conditions. Our goal is to use MetPy to calculate equivalent potential temperature ($\theta_e$), using the fields we have available in the RTMA.

This time we will also add some additional detail to our intermediate plot. To provide geographic reference, we will plot with a geographic projection and add coastlines.

Step 3a: Obtain Model Output ¶

Top

Step 3b: Prepare Model Output ¶

Top

Before we calculate any derived variables, we need to look up the required input variables. Looking at the metpy.calc documentation equivalent potential temperature, we can see that this function requires three variables: surface pressure, temperature, and dewpoint. Each of these variables are contained in our rtma_data Dataset, but we need to know what variable name each are stored under.

Now we pull out the fields (variables) we need for the $\theta_e$ calculation: surface pressure, surface temperature, and surface dewpoint. Here we are making use of xarray and MetPy to select the particular time we want.

Notice the use of the squeeze() method. This helps eliminate any stray dimensions for the vertical since we want a 2D array as output (some values are at a particular height, like 2 m, rather than the surface).

Now calculate equivalent potential temperature using the metpy.calc documentation

Recall using cartopy's crs (cartographic reference systems) projections from Day 1 of this workshop.

We can use MetPy shortcuts to create one of these objects for our data. This will help us during plotting, and can be used to change the geographic projection of our data on demand. This is the final preparation step we will take before plotting our data.

Step 3c: Visualize Model Output ¶

Top

The data are now prepared for plotting. Let's now use a different plotting technique, filled contours, for visualizing the $\theta_e$ field. Let's also project the field into the Robinson projection.

TIP: If you are working with xarray dataarrays parsed via metpy.parse_cf, and you aren't sure what your horizontal variables (e.g. x or lon or latitude) are named, or you just want a shortcut, you can access these quickly with xarray.metpy.x and xarray.metpy.y

For example, these two lines of code are equivalent for dataarrays parsed with metpy.parse.cf:

ax.contourf(theta_e['x'], theta_e['y'], theta_e, transform=rtma_crs)

ax.contourf(theta_e.metpy.x, theta_e.metpy.y, theta_e, transform=rtma_crs)

Step 4: Surface Observations ¶

Top

Finally, we will add station plots to our map. Creating a map of station plots is a bit different from our previous two map layers. We're no longer looking at 2D arrays of data, and instead we're looking at a collection of points containing several observations. Despite those differences, though, the same general workflow applies: obtain, prepare, and plot the data. However, the details of each will differ greatly from the satellite and $\theta_e$ plots.

Step 4a: Obtain Surface Observations ¶

Top

Step 4b: Prepare Surface Observations ¶

Top

metar_text now holds an object containing our METAR data in ascii text. Now we will prepare the data for plotting by converting it into a pandas DataFrame object using parse_metar_file(). A pandas DataFrame is similar to a spreadsheet in that the data are stored into rows and columns.

Example:

This structure allows us to more easily apply some quality control. Using Pandas, we will remove stations that do not have a useful lat/lon value (empty or null), as well as set missing current weather observations to an empty string.

Now we convert the Pandas DataFrame object into a dictionary object, where each column is now a unique key. This structure is what MetPy requires for creating the final station plot.

The next preparation step is to prepare the current weather conditions. METARs use text abbreviations to denote the current weather conditions, but MetPy requires that these codes be converted to the World Meteorological Organization (WMO) numerical code for plotting.

For example:

We will do this conversion by passing the current_wx1 key from the sfc_data dictionary through the wx_code_to_numeric() MetPy function.

Step 4c: Visualize Surface Observations ¶

Top

Now that we have completed all of our data preparation, we can start building the station plot. To do this, we need MetPy's StationPlot class. The StationPlot class acts as a template that we fill in with the surface observations stored in the sfc_data variable. The station plot is broken up into regions surrounding a center point. The variables we'll be plotting are:

1. air temperature (red, north west NW region)
2. dew point temperature (blue, south west SW region)
3. cloud coverage (white, center C region)
4. current weather symbol/wx code (blue, east E region)
5. wind speed and direction (wind barb).

Example:

NOTE: We are making a few changes from the standard station plot for accessibility and readability. Traditionally, the dewpoint value is displayed in green. However, we are making a modification to make our final plot more accessible to colorblindness. We are also moving the wx code to the east (as oppposed to west) to increase readability.

Oh. Oh no.

Because there are so many surface observations available, we ended up with an illegible plot. Let's now select only a subset of the stations to plot by creating a mask (filter) over our surface observations. The output should limit the number of stations to one station every 175 km.

Step 5: Create multi-layer plot ¶

Top

Now that we have obtained, prepared, and visualized all of the individual datasets, we can now put them all together in a multi-layer plot. Many of the visualization tools you previously used will be again applied here, with the addition of a few other tips to help with the readability of the final map.