This is a basic introduction to real-world data analysis and visualization with Python. To see the finished product and a bit more analysis, read my blog.
The first step in data analysis is finding the right data. For this project we'll be using the Global Historical Climatology Network (GHCN) monthly data compiled by the US government. In a previous post I explain where this data comes from.
A professor of mine, and an experienced data scientist, once mentioned that "80% of the effort in machine learning is getting the source data into a manageable format." The following section simply does exactly that. :) The source files are in a custom fixed-width-column format specified on the NOAA website.
The first block here parses three *.dat files, one each for minimum, maximum and average temperature measurements. Each of these files consists of one line per station per year. A line includes the station ID, year of observation, and measurements for all 12 months.
These files are simply read in, split at the correct locations, and written out to a CSV file. This will allow us to import the data into pandas directly using pandas's read_csv functionality.
In [1]:
from __future__ import division
from glob import glob
import csv,re
for name in glob('data/raw/*.dat'):
label = name[15:19]
print "Parsing data:",label
with open(name) as in_file:
with open("data/csv/"+label+".csv",'w') as out_file:
writer = csv.writer(out_file)
header = ["id","year"]
header += [str(month+1)+"_"+label for month in xrange(12)]
writer.writerow(header)
for line in in_file:
# if int(line[11:15]) < 2000:
# continue
components = [line[0:11],line[11:15]]
components += [int(line[19+x*8:24+x*8])/100 for x in range(12)]
writer.writerow(components)
A country-codes file is provided, which correlates each station ID with a country. Because this file is small we can read it straight into a python dictionary of the form {[country_code]: [country_name]}, which we will use in the next step to associate countries with station IDs.
In [2]:
def title(name):
"""Title case with correct handling of internal apostrophes"""
return re.sub("([a-z])'([A-Z])", lambda m: m.group(0).lower(), name.title())
cc_file = open('data/raw/country-codes')
ccodes = {l[:3]:title(l[4:].strip()) for l in cc_file}
Index files, ending in *.inv, correlate the station IDs in the .dat files with station metadata, such as latitude, longitude and station name. The following lines open each .inv file and convert them to a CSV, as above. It also uses the ccodes index from above to look up the country name for each file.
In [3]:
for name in glob('data/raw/*.inv'):
label = name[15:19]
print "Parsing index:",label
with open(name) as in_file:
with open("data/csv/"+label+"-inv.csv",'w') as out_file:
writer = csv.writer(out_file)
writer.writerow(["id","lat","long","name","country"])
for line in in_file:
# if not line[74:79].strip().isdigit() or int(line[74:79]) < 100:
# continue
name = title(line[38:68])
name = name[:name.find(" ")].strip()
components = [line[0:11],line[12:20],line[21:30],name,ccodes[line[:3]]]
writer.writerow(components)
Now that all the relevant data is stored as easily-parsed CSV files, we can begin processing it!
The main tool we'll be using for this is a pandas DataFrame. Think of a dataframe as much like a table in SQL or even a spreadsheet in Excel with lots of data indexing, retreival and aggregation tools built right in. For more information on DataFrames, see the pandas introduction. Pandas lets you build a dataframe from a dictionary, but in this case we'll be reading one right from a CSV file on disk. It looks something like this:
In [4]:
import pandas as pd
tavg_dataframe = pd.read_csv("data/csv/tavg.csv")
tavg_dataframe.head()
Out[4]:
Now that we're able to open our CSV files and create dataframes, we can do some basic number crunching with them. For this project we're not interested in the month-by-month breakdown of temperatures, we just want the average for a year. So for each of our dataframes, we need to take the mean of the months. That's done with the line
typ_data[typ] = typ_data.drop(['id','year'],axis=1).mean(axis=1,skipna=True)
typ_data is a dataframe of the schema we can see above. For each row, we want to average together all the columns except year and id. That's why we call drop -- that removes those two columns for the rest of the calculation. We then take the mean of the rest of the columns, and store it back in the dataframe as typ (a stand-in for tavg, tmin or tmax). The argument axis=1 means that we want to operate across columns, instead of rows which is the default.
We then remove all columns except those we're interested in, with the syntax typ_data = typ_data[['id','year',typ]]. This assigns a new dataframe to typ_data composed only of the columns from the original included in the double brackets.
Next, we merge all three dataframes together one by one. This works just like a join in SQL, and gives you similar options. We need to make sure that all three readings are from the same year are matched together, so we join on id and year. By default, this performs an inner join, so only stations and years with all three types of readings are included in the final dataframe.
Finally, we merge the completed dataframe with the metadata index, giving us a latitude, longitude, station name and country for each data reading.
In [5]:
import numpy as np
data = None
for typ in ["tavg","tmin","tmax"]:
# -99.99 is used in this dataset to represent missing data, so we replace it with NaN (not a number)
# This keeps it from being used in calculations.
typ_data = pd.read_csv("data/csv/"+typ+".csv").replace(-99.99,np.NaN)
typ_data[typ] = typ_data.drop(['id','year'],axis=1).mean(axis=1,skipna=True)
typ_data = typ_data[['id','year',typ]]
if data is None:
data = typ_data
else:
data = pd.merge(data, typ_data, on=['year','id'])
index = pd.read_csv("data/csv/tavg-inv.csv")
data = pd.merge(data, index, on=['id'])
data = data.dropna()
data.head()
Out[5]:
Just for fun, let's take a quick look at the stations we've got.
This plot shows the number of stations with reported readings by year. That big drop in the 90's is interesting!
To produce this, we just need to combine our dataframe with some simple plotting code from matplotlib. All we need for a simple plot like this are two lists: "x" and "y" values. For the "x" values, we simply enumerate all of the years with available data. This is done with data['year'].min() and data['year'].max() -- the brackets in this case select the year column, and .min() and .max() perform the expected aggregation functions.
We then count the number of data points for each year, by calling data[data.year == year].count(). This selects all of the rows whose year matches the year in question, and then counts them. A simple query. :)
The matplotlib code is fairly self-explanatory -- we just pass in the two lists, label the axes, and show the plot.
In [6]:
import matplotlib.pyplot as plt
%matplotlib inline
years = xrange(data['year'].min(),data['year'].max())
counts = [data[data.year == year].count() for year in years]
plt.plot(years,counts,color='red')
plt.xlabel('Year'); plt.ylabel('Stations Online')
plt.show()
The next section finds the average temperature for each station for the periods of 1950-1954 and 2008-2012. Most of these operations we've seen before, but we do introduce two new ones -- queries with multiple selectors and groupby.
early_years = data[(data.year >= 1950) & (data.year < 1955)][['id','tavg']]
This line uses simple ampersand syntax (&) to show that we want to make two selections on the data -- the year should be greater than or equal to 1950, and less than 1955. We also select only the id and tavg columns from the result.
early_years = early_years.groupby('id').mean()
We then use the powerful groupby syntax to find the mean of all the years in the collection that have the same station ID. This gives us the average temperature over the whole range for each station.
Finally, we use a simple subtraction operator to find the difference in averages. dropna is called to remove stations that do not have data available in both ranges.
In [7]:
from datetime import datetime
data = data[(data.year >= 1950) & (data.year <= 2013)]
data['date'] = data['year'].apply(lambda x: datetime(x,1,1))
early_years = data[(data.year >= 1950) & (data.year < 1955)][['id','tavg']]
early_years = early_years.groupby('id').mean()
late_years = data[(data.year >= 2008) & (data.year < 2013)][['id','tavg']]
late_years = late_years.groupby('id').mean()
deltas = (late_years['tavg']-early_years['tavg']).dropna()
We combine the evaluated deltas with the index using the same square brackets to assign a column.
The last line groups by country, and counts the occurances of each. 1219 of the 1529 stations represented are in the USA.
In [8]:
change = index.copy().set_index('id')
change['delta'] = deltas
change = change.dropna()
print "Total stations represented:",len(change)
change.groupby('country').count().sort('country',ascending=False)[['country']]
Out[8]:
describe gives you standard dataset summary data.
In [9]:
print change['delta'].describe()
matplotlib is great for quick analysis and visualizations. Let's plot a histogram of our temperature deltas.
In [10]:
plt.hist(list(change['delta']),50)
plt.xlim([-5,5])
plt.title('Average Temperature Change, 1950-2010')
plt.xlabel('Delta Temperature (C)')
plt.ylabel('Count of Stations')
plt.savefig('viz/temp_change.png')
plt.show()
The final step is mapping the data. There are three major steps to this process: (1) Setting the color of each marker, (2) creating the temperature over time graph for each marker, and (3) plotting the markers on the map.
We want to assign a color, blue or red, to each marker depending on its change over the period. The first step is to discretize the deltas into 7 equally-sized categories. Because most stations had a change of less than +/-5° C, we will use that as the limits of our scale. The first two lines find any changes of greater than 5° and cap them. After that, the cut method is used to discretize the delta column into 7 categories. I chose a 7-tone blue-to-red scale from ColorBrewer to represent the colors of the markers, which will be assigned based on the discretized delta.
In [11]:
change.delta[change.delta > 5] = 5
change.delta[change.delta < -5] = -5
change['cat'] = pd.cut(change['delta'], 7)
cats = change.sort('delta')['cat'].unique()
colors = ["#2166ac","#67a9cf","#d1e5f0","#f7f7f7","#fddbc7","#ef8a62","#b2182b"]
cat_colors = {cat:color for cat, color in zip(cats,colors)}
list(cats)
Out[11]:
Visualization is done with folium for the map and vincent for charts. These both have simple and well-documented APIs.
To create the graph for each station, we first select the readings that pertain to the station of interest, using the familiar bracket syntax for selection. The next line does several things.
weather_frame = weather_frame.rename(columns={'tmax':'Max','tmin':'Min','tavg':'Average'}).set_index('date')
We first have to rename the columns to human-friendly names like Max, Min and Average. Vincent will use these column names in its legend. Next, we set the dataframe index to date. This lets vincent know that we want a time-series graph with the date on the X axis.
Vincent doesn't need any more configuration than that, so we can simply create a line graph with chart = vincent.Line(weather_frame). The rest of the code performs fairly standard formatting of the chart, and is all found in the vincent API linked above.
In [12]:
import folium, vincent
vincent.core.initialize_notebook()
def popup_chart(station_id):
weather_frame = data[data.id == station_id][['date','tmax','tavg','tmin']]
weather_frame = weather_frame.rename(columns={'tmax':'Max','tmin':'Min','tavg':'Average'}).set_index('date')
station = index[index.id == station_id].iloc[0] # Find the station in the index
chart = vincent.Line(weather_frame)
chart.axis_titles(x=station['name']+", "+station["country"],y="Temp (C)")
chart.legend(title="Key")
chart.width=400
chart.height=200
chart.title=station['name']+", "+station["country"]
chart.padding={'top':20,'right':80,'bottom':40,'left':40}
chart.to_json("viz/historical_map/data_{0}.json".format(station['id']))
return chart
chart = popup_chart(42572734000)
chart.display()
With folium, creating the map itself is fairly simple! The configuration options, such as the default zoom level and tile set to use, are found in the folium docs.
We then iterate through all our stations, and add a marker for each one. The popup is set by generating a vincent chart in the above manner for each data point. The fill color is set with the previously generated cat_colors dictionary. We save the final map to disk, and we are good to go!
(Note: the final map won't show up in the online ipython notebook viewer because it relies on Javascript. However, you can see it here.)
In [13]:
folium.initialize_notebook()
f_map = folium.Map(
location=[0,0],
zoom_start=2,
tiles="Mapbox Bright")
for station_id,station in change.iterrows():
f_map.polygon_marker([station['lat'],station['long']],
radius=5,
popup=(popup_chart(station_id),"data_{0}.json".format(station_id)),
line_color='grey',
line_weight=1,
fill_color=cat_colors[station['cat']],
fill_opacity=1
)
f_map.create_map(path='viz/historical_map/map.html')
f_map
Out[13]: