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.

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