Atmospheric pressure in De Bilt

This entry is part 9 of 10 in the series NumPy Weather

Atmospheric pressure is the pressure exerted by air in the atmosphere. It is defined as force divided by area. The KNMI De Bilt data file has measurements in .1 hPa for the average, minimum and maximum daily pressures. We will plot a histogram of the average pressure and monthly minimums, maximums and averages.

Imports

We will import the NumPy (line 1) module, masked arrays NumPy module (line 2) and Matplotlib (line 3).

1
2
3
4
5
6
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
import sys
from datetime import datetime as dt
import calendar as cal

Loading the Data

We will load (line 3) the dates converted to months (line 2), average pressure, minimum and maximum pressure into NumPy arrays. Again missing values needed to be converted (line 1) into NaNs (not a number).

1
2
3
to_float = lambda x: 0.1 * float(x.strip() or np.nan)
to_month = lambda x: dt.strptime(x, "%Y%m%d").month
months, avg_p, max_p, min_p = np.loadtxt(sys.argv[1], delimiter=',', usecols=(1, 25, 26, 28), unpack=True, converters={1: to_month, 25: to_float, 26: to_float, 28: to_float})

Statistics

Values are missing from the pressure value columns, so we have to create masked arrays out of the NumPy arrays. The snippet below prints some simple statistics.

1
2
3
4
5
6
7
8
max_p = ma.masked_invalid(max_p)
print "Maximum Pressure", max_p.max()
 
avg_p = ma.masked_invalid(avg_p)
print "Average Pressure", avg_p.mean(), "Std Dev", avg_p.std()
 
min_p = ma.masked_invalid(min_p)
print "Minimum Pressure", min_p.max()

Which prints the following values.

Maximum Pressure 1050.4
Average Pressure 1015.14058231 Std Dev 9.85889134337
Minimum Pressure 1045.1

Monthly Aggregates

I compute monthly averages, minimums and maximums with the code below.

1
2
3
4
5
6
7
8
9
10
monthly_pressure = []
maxes = []
mins = []
month_range = np.arange(int(months.min()), int(months.max()))
 
for month in month_range:
   indices = np.where(month == months)
   monthly_pressure.append(avg_p[indices].mean())
   maxes.append(max_p[indices].max())
   mins.append(min_p[indices].min())


Plotting

We will draw a histogram (line 3) of the average daily pressures and associated Gaussian curve (line 6). In addition we will plot monthly aggregate values as prepared in the previous section.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
plt.subplot(211)
plt.title("Pressure Histogram")
a, bins, b = plt.hist(avg_p.compressed(), 200, normed=True)
stdev = avg_p.std()
avg = avg_p.mean()
plt.plot(bins, 1/(stdev * np.sqrt(2 * np.pi)) * np.exp(- (bins - avg)**2/(2 * stdev**2)), 'r-')
 
ax = plt.subplot(212)
plt.title("Monthly Pressure")
plt.plot(month_range, monthly_pressure, 'bo', label="Average")
plt.plot(month_range, maxes, 'r^', label="Maximum Values")
plt.plot(month_range, mins, 'g>', label="Minumum Values")
ax.set_xticklabels(cal.month_abbr[::2])
plt.legend(prop={'size':'x-small'}, loc='best')
ax.set_ylabel('hPa')
plt.show()

The following plots are produced.

 photo pressure_zpsc2bab245.png

As you can see the bell curve fits the distribution of average daily pressures almost perfectly. The monthly average pressure seems to be constant. I find this suspicious. I think that it must be an error. Today’s code listing is given below.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
import sys
from datetime import datetime as dt
import calendar as cal
 
 
to_float = lambda x: 0.1 * float(x.strip() or np.nan)
to_month = lambda x: dt.strptime(x, "%Y%m%d").month
months, avg_p, max_p, min_p = np.loadtxt(sys.argv[1], delimiter=',', usecols=(1, 25, 26, 28), unpack=True, converters={1: to_month, 25: to_float, 26: to_float, 28: to_float})
 
max_p = ma.masked_invalid(max_p)
print "Maximum Pressure", max_p.max()
 
avg_p = ma.masked_invalid(avg_p)
print "Average Pressure", avg_p.mean(), "Std Dev", avg_p.std()
 
min_p = ma.masked_invalid(min_p)
print "Minimum Pressure", min_p.max()
 
monthly_pressure = []
maxes = []
mins = []
month_range = np.arange(int(months.min()), int(months.max()))
 
for month in month_range:
   indices = np.where(month == months)
   monthly_pressure.append(avg_p[indices].mean())
   maxes.append(max_p[indices].max())
   mins.append(min_p[indices].min())
 
plt.subplot(211)
plt.title("Pressure Histogram")
a, bins, b = plt.hist(avg_p.compressed(), 200, normed=True)
stdev = avg_p.std()
avg = avg_p.mean()
plt.plot(bins, 1/(stdev * np.sqrt(2 * np.pi)) * np.exp(- (bins - avg)**2/(2 * stdev**2)), 'r-')
 
ax = plt.subplot(212)
plt.title("Monthly Pressure")
plt.plot(month_range, monthly_pressure, 'bo', label="Average")
plt.plot(month_range, maxes, 'r^', label="Maximum Values")
plt.plot(month_range, mins, 'g>', label="Minumum Values")
ax.set_xticklabels(cal.month_abbr[::2])
plt.legend(prop={'size':'x-small'}, loc='best')
ax.set_ylabel('hPa')
plt.show()

List for April 10, 2013

http://storify.com/inningPalmer/list-for-april-10-2013

Series NavigationMonthly Precipitation in De BiltAtmospheric Humidity in De Bilt
By the author of NumPy Beginner's Guide, NumPy Cookbook and Instant Pygame. If you enjoyed this post, please consider leaving a comment or subscribing to the RSS feed to have future articles delivered to your feed reader.
Share
This entry was posted in programming and tagged , . Bookmark the permalink.