Global Warming

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

According to the Global Warming theory the temperature on Earth has increased on average since the end of the nineteenth century. From the previous century until now the temperature supposedly has gained about 0.8 degrees. Apparently most of this warming has happened in the last two or three decades. In the future we can expect the temperature to rise even more leading to droughts, heat waves and other unpleasant phenomena. Obviously some regions will be harder hit than others. Several solutions have been proposed including reduction of greenhouse gas emissions and geoengineering by spreading special gases in the atmosphere in order to reflect more sunlight.

The following factors can contribute to changes in the Earth’s temperature:

  • greenhouse gases
  • the Sun
  • volcanic eruptions

 

Some Checks

The data I downloaded from the Dutch meteorological institute KNMI (see the previous posts in this series) is not sufficient to prove whether global warming is real or not, but we can certainly check some things. For instance, we can check whether the temperature in De Bilt ( that’s where the data was collected) in the first half of the data set is lower than in the second half. Another thing we can do is plot yearly average temperatures. De Bilt by the way, as far as I know, is a small town in Central Netherlands without heavy industry.

Imports

We need to import NumPy (line 1) and Matplotlib (line 2) to create plots later on.

1
2
3
4
import numpy as np
import matplotlib.pyplot as plt
import sys
from datetime import datetime as dt


Loading the data

We will load (line 3) the average daily temperatures and the corresponding dates. Actually we will convert the dates to years (line 1) immediately to be able to calculate yearly average temperatures.

1
2
3
4
5
6
7
8
9
10
11
to_year = lambda x: dt.strptime(x, "%Y%m%d").year
 
years, avg_temp = np.loadtxt(sys.argv[1], delimiter=',', usecols=(1, 11), unpack=True, converters={1: to_year})
 
# Measurements are in .1 degrees Celcius
avg_temp = .1 * avg_temp
 
N = len(avg_temp)
print "First Year", years[0], "Last Year", years[-1]
assert N == len(years)
assert years[:N/2].mean() < years[N/2:].mean()

As you can see some sanity checking occurs at the end of the snippet. Line 9 prints the following:

First Year 1901.0 Last Year 2013.0

The two halves

After dividing the average daily temperature values in two halves, we can calculate and compare the arithmetic means for both halves.

1
2
print "First half average", avg_temp[:N/2].mean(), "Std Dev", avg_temp[:N/2].std()
print "Second half average", avg_temp[N/2:].mean(), "Std Dev", avg_temp[N/2:].std()

This gives us:

First half average 9.19078446678 Std Dev 6.42457006016
Second half average 9.78066152795 Std Dev 6.34152195332

It seems that the average temperature is slightly higher in the second half of the dataset.

Yearly Average Temperatures
Computing the yearly average temperatures is trivial. The steps required are:

  1. Create a range for the years in the data set (line 2). We will omit the last year 2013, because we don’t have a full year’s worth of data for that year.
  2. For each year find the array indices corresponding to that year (line 5).
  3. Compute the mean for each year and store it (line 6).
1
2
3
4
5
6
7
8
9
10
11
avgs = []
year_range = range(int(years[0]), int(years[-1]) - 1)
 
for year in year_range:
   indices = np.where(years == year)
   avgs.append(avg_temp[indices].mean())
 
plt.plot(year_range, avgs, 'r-', label="Yearly Averages")
plt.plot(year_range, np.ones(len(avgs)) * np.mean(avgs))
plt.legend(prop={'size':'x-small'})
plt.show()

We get the plot below as a result. For comparison an average of all average temperatures is also drawn through the middle of the plot. Notice how the yearly average temperatures seem to be on the rise from 1980 on.

 photo globalWarming_zpsf3144167.png

Below is the code in its entirety.

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
import numpy as np
import matplotlib.pyplot as plt
import sys
from datetime import datetime as dt
 
to_year = lambda x: dt.strptime(x, "%Y%m%d").year
 
years, avg_temp = np.loadtxt(sys.argv[1], delimiter=',', usecols=(1, 11), unpack=True, converters={1: to_year})
 
# Measurements are in .1 degrees Celcius
avg_temp = .1 * avg_temp
 
N = len(avg_temp)
print "First Year", years[0], "Last Year", years[-1]
assert N == len(years)
assert years[:N/2].mean() < years[N/2:].mean()
print "First half average", avg_temp[:N/2].mean(), "Std Dev", avg_temp[:N/2].std()
print "Second half average", avg_temp[N/2:].mean(), "Std Dev", avg_temp[N/2:].std()
 
avgs = []
year_range = range(int(years[0]), int(years[-1]) - 1)
 
for year in year_range:
   indices = np.where(years == year)
   avgs.append(avg_temp[indices].mean())
 
plt.plot(year_range, avgs, 'r-', label="Yearly Averages")
plt.plot(year_range, np.ones(len(avgs)) * np.mean(avgs))
plt.legend(prop={'size':'x-small'})
plt.show()

Unfortunately this is all I have time for today. I am still doing research on meteorology, statistics and what not.


Archived on April 3, 2013

http://storify.com/inningPalmer/archived-on-april-3-2013

Series NavigationDaily Temperature RangeSun Radiation versus Temperature
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.