This is a more in-depth example based off of the basic meteogram, but using real data from the Oklahoma Mesonet. This data is from the El Reno, Oklahoma mesonet station on May 24, 2011, when the outer circulation of an EF-5 tornado passed directly over the mesonet station. This uses some advanced techniques, such as reading data from a file, combining multiple plots, defining functions, and more. While not all of these topics have been covered directly, this is to give you a more complete example of how a real meteogram would be generated.
As always, we will start by importing the libraries we need for this demo.
In [1]:
## we need to import a few libraries to get us going
import numpy as np
import matplotlib.pyplot as plt
In [2]:
## we will use this function to convert Relative Humidity to Dewpoint Temperature
def TD(RH, T):
'''
An approximation used to convert values of Temperature (C)
and Relative Humidity (%) into values of Dewpoint Temperature (C).
Source: http://andrew.rsmas.miami.edu/bmcnoldy/Humidity.html
'''
return 243.04*(np.log(RH/100)+((17.625*T)/(243.04+T)))/(17.625-np.log(RH/100)-((17.625*T)/(243.04+T)))
In this next cell, we will read in the data. For convenience, I am using a function from the numpy library that makes reading in data rather simple. However, in a future example, we will learn how to read in this data from scratch.
In [3]:
## open and read the file
## STID: Station ID
## STNM: Station Number
## TIME: Minutes after base time (typically 0000 UTC)
## RELH: Relative Humidity (%)
## TAIR: Air Temperature (C)
## WSPD: 5-minute averaged wind speed (m/s) at 10m.
## WVEC: 5-minute averaged wind velocity (m/s) (speed and direction accounted for) at 10m.
## WDIR: 5-minute averaged wind direction (degrees) at 10m.
## WDSD: Standard deviation of wind direction (degrees) during the 5-minute interval.
## WSSD: Standard deviation of wind speed (m/s) during the 5-minute interval.
## WMAX: Highest 3-second wind speed (m/s) at 10m sample.
## RAIN: Liquid precipitation accumulation (mm) since 0000 UTC. Frozen precipitation cannot be recorded
## until it melts; therefore, precipitation from snow may not be recorded until several days
## after the snow event.
## PRES: 5-minute averaged atmospheric pressure (mb).
## SRAD: 5-minute averaged downwelling global solar radiation (Watts/m^2).
## TA9M: 5-minute averaged air temperature (C) at 9m
## WS2M: 5-minute averaged wind speed (m/s) at 2m
## The rest of the variables are all related to soil moisture. The complete description of variables
## can be found at http://www.mesonet.org/index.php/site/about/mdf_mts_files
filename = '../data/20110524elre.mts'
STID, STNM, TIME, RELH, TAIR, WSPD, WVEC, WDIR, WDSD, \
WSSD, WMAX, RAIN, PRES, SRAD, TA9M, WS2M, TS10, TB10, \
TS05, TB05, TS30, TR05, TR25, TR60 = np.genfromtxt(filename, dtype=float, skip_header=3,
unpack=True, usemask=True, missing_values="-999")
## Dewpoint Temperature (C)
TDEW = TD(RELH, TAIR)
This cell is where the plotting will take place. We don't take ourselves very seriously at HOOT, so we will be having a little fun and using the xkcd function built into matplotlib. For those of you not familiar with xkcd, stop what you are doing, and go educate yourself right now... http://xkcd.com
For the record, I created this tutorial entirely off of code I ripped from pieces of the matplotlib example gallery. For your reference, this example gallery can be found here: http://matplotlib.org/gallery.html Nearly all of this code is spliced from pieces of example code.
In [4]:
## make our plot goofy
plt.xkcd()
## set up the plotting figure
fig, (ax0, ax1, ax2, ax3) = plt.subplots(figsize=(12,9), nrows=4)
## we are copying the fourth axis so that we can plot multiple things on it
ax3_copy = ax3.twinx()
## we're going to use this information to plot the hours on the x axis
timevals = np.arange(0, 1400, 60)
hours = np.arange(0, 24, 1)
## plot the main title of the whole image
plt.suptitle("Mesonet Funtimes!")
## set the first figure's axes
ax0.set_xlim(0, 1400)
ax0.set_ylim(10, 30)
## set the tick interval and tick labels for the x axis
ax0.set_xticks(timevals)
ax0.set_xticklabels(hours)
## set the second figure's axes
ax1.set_xlim(0, 1400)
ax1.set_ylim(935, 960)
## set the tick interval and tick labels for the x axis
ax1.set_xticks(timevals)
ax1.set_xticklabels(hours)
## set the third figure's axes
ax2.set_xlim(0, 1400)
ax2.set_ylim(0, 70)
## set the tick interval and tick labels for the x axis
ax2.set_xticks(timevals)
ax2.set_xticklabels(hours)
## set the fourth figure's axes
ax3.set_xlim(0, 1400)
ax3.set_ylim(0, 30)
## set the tick interval and tick labels for the x axis
ax3.set_xticks(timevals)
ax3.set_xticklabels(hours)
## create second axis for the third figure.
## this will be used to plot wind direction
## on the same frame.
ax3_copy.set_xlim(0, 1400)
ax3_copy.set_ylim(0, 360)
ax3_copy.set_xticks(timevals)
ax3_copy.set_xticklabels(hours)
## set the tick interval and the labels with it
ax3_copy.set_yticks((0, 90, 180, 270, 360))
ax3_copy.set_yticklabels(('N', 'E', 'S', 'W', 'N'))
## plot the first figure
## we want to draw a color fill under the line. For the temperature
## trace, we only want to fill where TAIR > TDEW
ax0.fill_between(TIME, TAIR, TDEW, where=TAIR>TDEW, facecolor="#E90000", alpha=0.3)
ax0.fill_between(TIME, TDEW, facecolor="#009A0D", alpha=0.3)
## then we want to plot a thick line to show the trace
ax0.plot(TIME, TAIR, '#E90000', linewidth=2.0)
ax0.plot(TIME, TDEW, '#009A0D', linewidth=2.0)
## then set a label for the plot
ax0.set_ylabel( "Temperature & Dewpoint (c)")
## plot the second figure
ax1.fill_between(TIME, PRES, facecolor="#705400", alpha=0.3)
ax1.plot(TIME, PRES, color="#705400", linewidth=2.0)
ax1.set_ylabel("Pressure (mb)")
## plot the third figure
ax2.fill_between(TIME, WMAX, facecolor="#0008A1", alpha=0.3)
ax2.plot(TIME, WMAX, color='#0008A1', linewidth=2.0)
ax2.set_ylabel("Max 10m Wind (m/s)")
## plot the fourth figure
## do the color fills
ax3.fill_between(TIME, WS2M, facecolor="#0A8CD2", alpha=0.3)
## only fill where WSPD > WS2M
ax3.fill_between(TIME, WSPD, where=WSPD>WS2M, facecolor="#004468", alpha=0.3)
## draw thicker line traces on top
ax3.plot(TIME, WS2M, color="#0A8CD2", linewidth=2.0)
ax3.plot(TIME, WSPD, color="#004468", linewidth=2.0)
## plot the wind direction as a scatter plot.
ax3_copy.scatter(TIME, WDIR, s=30, facecolors='none', edgecolors='y')
## set the labels
ax3.set_ylabel("2m and 10m Wind (m/s)")
ax3_copy.set_ylabel("Wind Direction")
## Let's annotate the tornado
## find the index of the array where the pressure is lowest
idx = np.where(PRES == PRES.min())
## get the time, pres (x,y) value at the minimum
x = TIME[idx]
y = PRES[idx]
## we want our text to sit at 700 minutes, and 965 mb
x2, y2 = 700, 965
## draw the annotation
ax1.annotate('SUCH PRESSURE DROP', xy=(x,y), arrowprops=dict(arrowstyle='->'), xytext=(x2, y2) )
idx = np.where( WMAX == WMAX.max())[0]
x = TIME[idx]
y = WMAX[idx]
x2, y2 = 700, 30
ax2.annotate('VERY WIND', xy=(x,y), arrowprops=dict(arrowstyle='->'), xytext=(x2, y2))
idx = np.where( WS2M == WS2M.max())[0]
x = TIME[idx]
y = WMAX[idx]
x2, y2 = 1000, 35
ax3.annotate('WOW', xy=(x,y), arrowprops=dict(arrowstyle='->'), xytext=(x2, y2))
## this function will get rid of whitespace
plt.tight_layout()
## show the plot and then clear it once the user closes it
plt.show()
plt.close()
In [6]: