Last updated: June 29th 2016
Welcome to a demo of Python's data analysis package called Pandas. Our goal is to learn about Data Analysis and transformation using Pandas while exploring datasets used to analyze climate change.
The global goal of this demo is to provide the tools to be able to try and reproduce some of the analysis done in the IPCC global climate reports published in the last decade (see for example https://www.ipcc.ch/pdf/assessment-report/ar5/syr/SYR_AR5_FINAL_full.pdf).
We are first going to load a few public datasets containing information about global temperature, global and local sea level infomation, and global concentration of greenhouse gases like CO2, to see if there are correlations and how the trends are to evolve, assuming no fundamental change in the system. For all these datasets, we will download them, visualize them, clean them, search through them, merge them, resample them, transform them and summarize them.
In the process, we will learn about:
1. Loading data
2. Pandas datastructures
3. Cleaning and formatting data
4. Basic visualization
5. Accessing data
6. Working with dates and times
7. Transforming datasets
8. Statistical analysis
9. Data agregation and summarization
10. Correlations and regressions
11. Predictions from auto regression models
In [1]:
%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
pd.set_option("display.max_rows", 16)
LARGE_FIGSIZE = (12, 8)
In [2]:
# Change this cell to the demo location on YOUR machine
%cd ~/Projects/pandas_tutorial/climate_timeseries/
%ls
More details, see http://pandas.pydata.org/pandas-docs/stable/io.html
To find all reading functions in pandas, ask ipython's tab completion:
In [3]:
#pd.read_<TAB>
In [4]:
pd.read_table?
Let's first load some temperature data which covers all lattitudes. Since read_table is supposed to do the job for a text file, let's just try it:
In [5]:
filename = "data/temperatures/annual.land_ocean.90S.90N.df_1901-2000mean.dat"
full_globe_temp = pd.read_table(filename)
full_globe_temp
Out[5]:
There is only 1 column! Let's try again stating that values are separated by any number of spaces:
In [6]:
full_globe_temp = pd.read_table(filename, sep="\s+")
full_globe_temp
Out[6]:
There are columns but the column names are 1880 and -0.1591!
In [7]:
full_globe_temp = pd.read_table(filename, sep="\s+", names=["year", "mean temp"])
full_globe_temp
Out[7]:
Since we only have 2 columns, one of which would be nicer to access the data (the year of the record), let's try using the index_col option:
In [8]:
full_globe_temp = pd.read_table(filename, sep="\s+", names=["year", "mean temp"],
index_col=0)
full_globe_temp
Out[8]:
Last step: the index is made of dates. Let's make that explicit:
In [9]:
full_globe_temp = pd.read_table(filename, sep="\s+", names=["year", "mean temp"],
index_col=0, parse_dates=True)
full_globe_temp
Out[9]:
Since every dataset can contain mistakes, let's load a different file with temperature data. NASA's GISS dataset is written in chunks: look at it in data/temperatures/GLB.Ts+dSST.txt
In [10]:
giss_temp = pd.read_table("data/temperatures/GLB.Ts+dSST.txt", sep="\s+", skiprows=7,
skip_footer=11, engine="python")
giss_temp
Out[10]:
QUIZ: What happens if you remove the skiprows? skipfooter? engine?
EXERCISE: Load some readings of CO2 concentrations in the atmosphere from the data/greenhouse_gaz/co2_mm_global.txt data file.
In [11]:
# Your code here
So far, we have only loaded temperature datasets. Climate change also affects the sea levels on the globe. Let's load some datasets with the sea levels. The university of colorado posts updated timeseries for mean sea level globably, per hemisphere, or even per ocean, sea, ... Let's download the global one, and the ones for the northern and southern hemisphere.
That will also illustrate that to load text files that are online, there is no more work than replacing the filepath by a URL n read_table:
In [12]:
# Local backup: data/sea_levels/sl_nh.txt
northern_sea_level = pd.read_table("http://sealevel.colorado.edu/files/current/sl_nh.txt",
sep="\s+")
northern_sea_level
Out[12]:
In [13]:
# Local backup: data/sea_levels/sl_sh.txt
southern_sea_level = pd.read_table("http://sealevel.colorado.edu/files/current/sl_sh.txt",
sep="\s+")
southern_sea_level
Out[13]:
In [14]:
# The 2015 version of the global dataset:
# Local backup: data/sea_levels/sl_ns_global.txt
url = "http://sealevel.colorado.edu/files/2015_rel2/sl_ns_global.txt"
global_sea_level = pd.read_table(url, sep="\s+")
global_sea_level
Out[14]:
There are clearly lots of cleanup to be done on these datasets. See below...
To be able to grab more local data about mean sea levels, we can download and extract data about mean sea level stations around the world from the PSMSL (http://www.psmsl.org/). Again to download and parse all tables in a webpage, just give read_html the URL to parse:
In [15]:
# Needs `lxml`, `beautifulSoup4` and `html5lib` python packages
# Local backup in data/sea_levels/Obtaining Tide Gauge Data.html
table_list = pd.read_html("http://www.psmsl.org/data/obtaining/")
In [16]:
# there is 1 table on that page which contains metadata about the stations where
# sea levels are recorded
local_sea_level_stations = table_list[0]
local_sea_level_stations
Out[16]:
That table can be used to search for a station in a region of the world we choose, extract an ID for it and download the corresponding time series with the URL http://www.psmsl.org/data/obtaining/met.monthly.data/< ID >.metdata
For more details, see http://pandas.pydata.org/pandas-docs/stable/dsintro.html
Now that we have used read_** functions to load datasets, we need to understand better what kind of objects we got from them to learn to work with them.
In [17]:
# Type of the object?
type(giss_temp)
Out[17]:
In [18]:
# Internal nature of the object
print(giss_temp.shape)
print(giss_temp.dtypes)
Descriptors for the vertical axis (axis=0)
In [19]:
giss_temp.index
Out[19]:
Descriptors for the horizontal axis (axis=1)
In [20]:
giss_temp.columns
Out[20]:
A lot of information at once including memory usage:
In [21]:
giss_temp.info()
A series can be constructed with the pd.Series constructor (passing a list or array of values) or from a DataFrame, by extracting one of its columns.
In [22]:
# Do we already have a series for the full_globe_temp?
type(full_globe_temp)
Out[22]:
In [23]:
# Since there is only one column of values, we can make this a Series without
# loosing information:
full_globe_temp = full_globe_temp["mean temp"]
Core attributes/information:
In [24]:
print(type(full_globe_temp))
print(full_globe_temp.dtype)
print(full_globe_temp.shape)
print(full_globe_temp.nbytes)
Probably the most important attribute of a Series or DataFrame is its index since we will use that to, well, index into the structures to access te information:
In [25]:
full_globe_temp.index
Out[25]:
It is always possible to fall back to a good old NumPy array to pass on to scientific libraries that need them: SciPy, scikit-learn, ...
In [26]:
full_globe_temp.values
Out[26]:
In [27]:
type(full_globe_temp.values)
Out[27]:
DataFrames can also be created manually, by grouping several Series together. Let's make a new frame from the 3 sea level datasets we downloaded above. They will be displayed along the same index. Wait, does that makes sense to do that?
In [28]:
# Are they aligned?
southern_sea_level.year == northern_sea_level.year
Out[28]:
In [29]:
# So, are they aligned?
np.all(southern_sea_level.year == northern_sea_level.year)
Out[29]:
So the northern hemisphere and southern hemisphere datasets are aligned. What about the global one?
In [30]:
len(global_sea_level.year) == len(northern_sea_level.year)
Out[30]:
For now, let's just build a DataFrame with the 2 hemisphere datasets then. We will come back to add the global one later...
In [31]:
mean_sea_level = pd.DataFrame({"northern_hem": northern_sea_level["msl_ib(mm)"],
"southern_hem": southern_sea_level["msl_ib(mm)"],
"date": northern_sea_level.year})
mean_sea_level
Out[31]:
Note: there are other ways to create DataFrames manually, for example from a 2D numpy array.
There is still the date in a regular column and a numerical index that is not that meaningful. We can specify the index of a DataFrame at creation. Let's try:
In [32]:
mean_sea_level = pd.DataFrame({"northern_hem": northern_sea_level["msl_ib(mm)"],
"southern_hem": southern_sea_level["msl_ib(mm)"]},
index = northern_sea_level.year)
mean_sea_level
Out[32]:
Now the fact that it is failing show that Pandas does auto-alignment of values: for each value of the index, it searches for a value in each Series that maps the same value. Since these series have a dumb numerical index, no values are found.
Since we know that the order of the values match the index we chose, we can replace the Series by their values only at creation of the DataFrame:
In [33]:
mean_sea_level = pd.DataFrame({"northern_hem": northern_sea_level["msl_ib(mm)"].values,
"southern_hem": southern_sea_level["msl_ib(mm)"].values},
index = northern_sea_level.year)
mean_sea_level
Out[33]:
The datasets that we obtain straight from the reading functions are pretty raw. A lot of pre-processing can be done during data read but we haven't used all the power of the reading functions. Let's learn to do a lot of cleaning and formatting of the data.
The GISS temperature dataset has a lot of issues too: useless numerical index, redundant columns, useless rows, placeholder (****) for missing values, and wrong type for the columns. Let's fix all this:
In [34]:
# The columns of the local_sea_level_stations aren't clean: they contain spaces and dots.
local_sea_level_stations.columns
Out[34]:
In [35]:
# Let's clean them up a bit:
local_sea_level_stations.columns = [name.strip().replace(".", "")
for name in local_sea_level_stations.columns]
local_sea_level_stations.columns
Out[35]:
We can also rename an index by setting its name. For example, the index of the mean_sea_level dataFrame could be called date since it contains more than just the year:
In [36]:
mean_sea_level.index.name = "date"
mean_sea_level
Out[36]:
In the full globe dataset, -999.00 was used to indicate that there was no value for that year. Let's search for all these values and replace them with the missing value that Pandas understand: np.nan
In [37]:
full_globe_temp == -999.000
Out[37]:
In [38]:
full_globe_temp[full_globe_temp == -999.000] = np.nan
full_globe_temp.tail()
Out[38]:
In [39]:
# We didn't set a column number of the index of giss_temp, we can do that afterwards:
giss_temp = giss_temp.set_index("Year")
giss_temp.head()
Out[39]:
In [40]:
# 1 column is redundant with the index:
giss_temp.columns
Out[40]:
In [41]:
# Let's drop it:
giss_temp = giss_temp.drop("Year.1", axis=1)
giss_temp
Out[41]:
In [42]:
# We can also just select the columns we want to keep:
giss_temp = giss_temp[[u'Jan', u'Feb', u'Mar', u'Apr', u'May', u'Jun', u'Jul',
u'Aug', u'Sep', u'Oct', u'Nov', u'Dec']]
giss_temp
Out[42]:
In [43]:
# Let's remove all these extra column names (Year Jan ...). They all correspond to the index "Year"
giss_temp = giss_temp.drop("Year")
giss_temp
Out[43]:
Let's also set **** to a real missing value (np.nan). We can often do it using a boolean mask, but that may trigger pandas warning. Another way to assign based on a boolean condition is to use the where method:
In [44]:
#giss_temp[giss_temp == "****"] = np.nan
giss_temp = giss_temp.where(giss_temp != "****", np.nan)
In [45]:
giss_temp.tail()
Out[45]:
While building the mean_sea_level dataFrame earlier, we didn't include the values from global_sea_level since the years were not aligned. Adding a column to a dataframe is as easy as adding an entry to a dictionary. So let's try:
In [46]:
mean_sea_level["mean_global"] = global_sea_level["msl_ib_ns(mm)"]
mean_sea_level
Out[46]:
The column is full of NaNs again because the auto-alignment feature of Pandas is searching for the index values like 1992.9323 in the index of global_sea_level["msl_ib_ns(mm)"] series and not finding them. Let's set its index to these years so that that auto-alignment can work for us and figure out which values we have and not:
In [47]:
global_sea_level = global_sea_level.set_index("year")
global_sea_level["msl_ib_ns(mm)"]
Out[47]:
In [48]:
mean_sea_level["mean_global"] = global_sea_level["msl_ib_ns(mm)"]
mean_sea_level
Out[48]:
EXERCISE: Create a new series containing the average of the 2 hemispheres minus the global value to see if that is close to 0. Work inside the mean_sea_level dataframe first. Then try with the original Series to see what happens with data alignment while doing computations.
In [49]:
# Your code here
Now that the sea levels are looking pretty good, let's got back to the GISS temperature dataset. Because of the labels (strings) found in the middle of the timeseries, every column only assumed to contain strings (didn't convert them to floating point values):
In [50]:
giss_temp.dtypes
Out[50]:
That can be changed after the fact (and after the cleanup) with the astype method of a Series:
In [51]:
giss_temp["Jan"].astype("float32")
Out[51]:
In [52]:
for col in giss_temp.columns:
giss_temp.loc[:, col] = giss_temp[col].astype(np.float32)
An index has a dtype just like any Series and that can be changed after the fact too.
In [53]:
giss_temp.index.dtype
Out[53]:
For now, let's change it to an integer so that values can at least be compared properly. We will learn below to change it to a datetime object.
In [54]:
giss_temp.index = giss_temp.index.astype(np.int32)
Removing missing values - once they have been converted to np.nan - is very easy. Entries that contain missing values can be removed (dropped), or filled with many strategies.
In [55]:
full_globe_temp
Out[55]:
In [56]:
full_globe_temp.dropna()
Out[56]:
In [57]:
# This will remove any year that has a missing value. Use how='all' to keep partial years
giss_temp.dropna(how="any").tail()
Out[57]:
In [58]:
giss_temp.fillna(value=0).tail()
Out[58]:
In [59]:
# This fills them with the previous year. See also temp3.interpolate
giss_temp.fillna(method="ffill").tail()
Out[59]:
Let's also mention the .interpolate method on a Series:
In [60]:
giss_temp.Aug.interpolate().tail()
Out[60]:
For now, we will leave the missing values in all our datasets, because it wouldn't be meaningful to fill them.
EXERCISE: Go back to the reading functions, and learn more about other options that could have allowed us to fold some of these pre-processing steps into the data loading.
Now they have been formatted, visualizing your datasets is the next logical step and is trivial with Pandas. The first thing to try is to invoke the .plot to generate a basic visualization (uses matplotlib under the covers).
In [61]:
full_globe_temp.plot()
Out[61]:
In [62]:
giss_temp.plot(figsize=LARGE_FIGSIZE)
Out[62]:
In [63]:
mean_sea_level.plot(subplots=True, figsize=(16, 12));
In [64]:
# Distributions of mean sean level globally and per hemisphere?
mean_sea_level.plot(kind="kde", figsize=(12, 8))
Out[64]:
QUIZ: How to list the possible kinds of plots that the plot method can allow?
In [65]:
# Distributions of temperature in each month since 1880
giss_temp.boxplot();
There are more plot options inside pandas.tools.plotting:
In [66]:
# Is there correlations between the northern and southern sea level timeseries we loaded?
from pandas.tools.plotting import scatter_matrix
scatter_matrix(mean_sea_level, figsize=LARGE_FIGSIZE);
We will confirm the correlations we think we see further down...
EXERCISE: Refer to exercises/aapl_adj_close_plot/aapl_adj_close_plot.ipynb
The general philosophy for accessing values inside a Pandas datastructure is that, unlike a numpy array that only allows to index using integers a Series allows to index with the values inside the index. That makes the code more readable.
In [67]:
full_globe_temp
Out[67]:
In [68]:
# By default [] on a series accesses values using the index, not the location in the series
# print(temp1[0]) # This would to fail!!
In [69]:
# This index is non-trivial though (will talk more about these datetime objects further down):
full_globe_temp.index.dtype
Out[69]:
In [70]:
first_date = full_globe_temp.index[0]
first_date == pd.Timestamp('1880')
Out[70]:
In [71]:
# By default [] on a series accesses values using the index, not the location in the series
print(full_globe_temp[pd.Timestamp('1880')])
# print(temp1[0]) # This would fail!!
In [72]:
# Another more explicit way to do the same thing is to use loc
print(full_globe_temp.loc[pd.Timestamp('1990')])
print(full_globe_temp.iloc[0], full_globe_temp.iloc[-1])
In [73]:
# Year of the last record?
full_globe_temp.index[-1]
Out[73]:
In [74]:
# New records can be added:
full_globe_temp[pd.Timestamp('2011')] = np.nan
In [75]:
# In 2D, same idea, though in a DF [] accesses columns (Series)
giss_temp["Jan"]
Out[75]:
In [76]:
# while .loc and .iloc allow to access individual values, slices or masked selections:
print(giss_temp.loc[1979, "Dec"])
In [77]:
# Slicing can be done with .loc and .iloc
print(giss_temp.loc[1979, "Jan":"Jun"]) # Note that the end point is included unlike NumPy!!!
print(giss_temp.loc[1979, ::2])
In [78]:
# Masking can also be used in one or more dimensions. For example, another way to grab every other month for the first year:
mask = [True, False] * 6
print(giss_temp.iloc[0, mask])
print(giss_temp.loc[1880, mask])
In [79]:
# We could also add a new column like a new entry in a dictionary
giss_temp["totals"] = giss_temp.sum(axis=1)
giss_temp
Out[79]:
In [80]:
# Let's remove this new column, we will learn to do this differently
giss_temp = giss_temp.drop("totals", axis=1)
More complex queries rely on the same concepts. For example what are the names, and IDs of the sea level stations in the USA?
In [81]:
local_sea_level_stations.columns
Out[81]:
In [82]:
american_stations = local_sea_level_stations["Country"] == "USA"
local_sea_level_stations.loc[american_stations, ["ID", "Station Name"]]
Out[82]:
More details at http://pandas.pydata.org/pandas-docs/stable/timeseries.html
Let's work some more with full_globe_temp's index since we saw it is special.
In [83]:
# Its dtype is NumPy's new 'datetime64[ns]':
full_globe_temp.index.dtype
Out[83]:
The advantage of having a real datetime index is that operations can be done efficiently on it. Let's add a flag to signal if the value is before or after the great depression's black Friday:
In [84]:
black_friday = pd.to_datetime('1929-10-29')
full_globe_temp.index > black_friday
Out[84]:
In [85]:
# Convert its index from timestamp to period: it is more meaningfull since it was measured and averaged over the year...
full_globe_temp.index = full_globe_temp.index.to_period()
full_globe_temp
Out[85]:
See also to_timestamp to conver back to timestamps and its how method to specify when inside the range to set the timestamp.
Another thing that can be done is to resample the series, downsample or upsample. Let's see the series converted to 10 year blocks or upscale to a monthly series:
In [86]:
# Frequencies can be specified as strings: "us", "ms", "S", "T", "H", "D", "B", "W", "M", "A", "3min", "2h20", ...
# More aliases at http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases
full_globe_temp.resample("M").mean()
Out[86]:
In [87]:
full_globe_temp.resample("10A").mean()
Out[87]:
The index for giss_temp isn't an instance of datetimes so we may want to generate such DatetimeIndex objects. This can be done with date_range and period_range:
In [88]:
# Can specify a start date and a number of values desired. By default it will assume an interval of 1 day:
pd.date_range('1/1/1880', periods=4)
Out[88]:
In [89]:
# Can also specify a start and a stop date, as well as a frequency
pd.date_range('1/1/1880', '1/1/2016', freq="A")
Out[89]:
Note that "A" by default means the end of the year. Other times in the year can be specified with "AS" (start), "A-JAN" or "A-JUN". Even more options can be imported from pandas.tseries.offsets:
In [90]:
from pandas.tseries.offsets import YearBegin
pd.date_range('1/1/1880', '1/1/2015', freq=YearBegin())
Out[90]:
Actually we will convert that dataset to a 1D dataset, and build a monthly index, so lets build a monthly period index
In [91]:
giss_temp_index = pd.period_range('1/1/1880', '12/1/2015', freq="M")
giss_temp_index
Out[91]:
Let's look at our local_sea_level_stations dataset some more, to learn more about it and also do some formatting. What is the range of dates and lattitudes we have, the list of countries, the range of variations, ...
In [92]:
# What about the range of dates?
local_sea_level_stations["Date"].min(), local_sea_level_stations["Date"].max(), local_sea_level_stations["Date"].iloc[-1]
Out[92]:
In [93]:
local_sea_level_stations.dtypes
Out[93]:
We don't see the range of dates because the dates are of dtype "Object", (usually meaning strings). Let's convert that using apply:
In [94]:
local_sea_level_stations["Date"].apply(pd.to_datetime)
Out[94]:
This apply method is very powerful and general. We have used it to do something we could have done with astype, but any custom function can be provided to apply.
In [95]:
local_sea_level_stations["Date"] = local_sea_level_stations["Date"].apply(pd.to_datetime)
# Now we can really compare the dates, and therefore get a real range:
print(local_sea_level_stations["Date"].min(), local_sea_level_stations["Date"].max())
EXERCISE: Use the apply method to search through the stations names for a station in New York. What is the ID of the station?
In [96]:
# Your code here
Now that we know the range of dates, to look at the data, sorting it following the dates is done with sort:
In [97]:
local_sea_level_stations.sort_values(by="Date")
Out[97]:
Since many stations last updated on the same dates, it is logical to want to sort further, for example, by Country at constant date:
In [98]:
local_sea_level_stations.sort_values(by=["Date", "Country"], ascending=False)
Out[98]:
Let's look at the GISS dataset differently. Instead of seeing the months along the axis 1, and the years along the axis 0, it would could be good to convert these into an outer and an inner axis along only 1 time dimension.
Stacking and unstacking allows to convert a dataframe into a series and vice-versa:
In [99]:
giss_temp.unstack?
unstacked = giss_temp.unstack()
unstacked
Out[99]:
In [100]:
# Note the nature of the result:
type(unstacked)
Out[100]:
The result is grouped in the wrong order since it sorts first the axis that was unstacked. Another transformation that would help us is transposing...
In [101]:
giss_temp.transpose()
Out[101]:
In [102]:
giss_temp_series = giss_temp.transpose().unstack()
giss_temp_series.name = "Temp anomaly"
giss_temp_series
Out[102]:
In [103]:
# Note the nature of the resulting index:
giss_temp_series.index
Out[103]:
In [104]:
# It is an index made of 2 columns. Let's fix the fact that one of them doesn't have a name:
giss_temp_series.index = giss_temp_series.index.set_names(["year", "month"])
In [105]:
# We can now access deviations by specifying the year and month:
giss_temp_series[1980, "Jan"]
Out[105]:
But this new multi-index isn't very good, because is it not viewed as 1 date, just as a tuple of values:
In [106]:
giss_temp_series.plot(figsize=LARGE_FIGSIZE)
Out[106]:
To improve on this, let's reuse an index we generated above with date_range:
In [107]:
giss_temp_series.index = giss_temp_index
giss_temp_series.plot(figsize=LARGE_FIGSIZE)
Out[107]:
Let's go back to the dataframe version of the GISS temperature dataset temporarily to analyze anomalies month per month. Like most functions on a dataframe, stats functions are computed column per column. They also ignore missing values:
In [108]:
monthly_averages = giss_temp.mean()
monthly_averages
Out[108]:
It is possible to apply stats functions across rows instead of columns using the axis keyword (just like in NumPy).
In [109]:
yearly_averages = giss_temp.mean(axis=1)
yearly_averages
Out[109]:
describe provides many descriptive stats computed at once:
In [110]:
mean_sea_level.describe()
Out[110]:
Let's remove high frequency signal and extract the trend:
In [111]:
full_globe_temp.plot()
rolled_series = full_globe_temp.rolling(window=10, center=False)
print(rolled_series)
rolled_series.mean().plot(figsize=LARGE_FIGSIZE)
Out[111]:
In [112]:
# To see what all can be done while rolling,
#pd.rolling_<TAB>
Let's look at our local_sea_level_stations dataset some more:
In [113]:
local_sea_level_stations.describe()
Out[113]:
.describe() only displays information about continuous Series. What about categorical ones?
In [114]:
local_sea_level_stations.columns
Out[114]:
In [115]:
local_sea_level_stations["Country"]
Out[115]:
In [116]:
local_sea_level_stations["Country"].describe()
Out[116]:
In [117]:
# List of unique values:
local_sea_level_stations["Country"].unique()
Out[117]:
In [118]:
local_sea_level_stations["Country"].value_counts()
Out[118]:
In [119]:
# To save memory, we can convert it to a categorical column:
local_sea_level_stations["Country"] = local_sea_level_stations["Country"].astype("category")
We can also create categorical series from continuous ones with the cut function:
In [120]:
categorized = pd.cut(full_globe_temp, 3, labels=["L", "M", "H"])
categorized
Out[120]:
In [121]:
# The advantage is that we can use labels and control the order they should be treated in (L < M < H)
categorized.cat.categories
Out[121]:
QUIZ: How much memory did we save? What if it was categorized but with dtype object instead of category?
Now that we have a good grasp on our datasets, Let's transform and analyze them some more to prepare them to compare them. The 2 function(alities)s to learn about here are groupby and pivot_table.
Let's explore the sea levels, first splitting into calendar years to compute average sea levels for each year:
In [122]:
mean_sea_level
Out[122]:
In [123]:
mean_sea_level = mean_sea_level.reset_index()
mean_sea_level
Out[123]:
In [124]:
# Groupby with pandas can be done on a column or by applying a custom function to the index.
# If we want to group the data by year, we can build a year column into the DF:
mean_sea_level["year"] = mean_sea_level["date"].apply(int)
mean_sea_level
Out[124]:
In [125]:
sl_grouped_year = mean_sea_level.groupby("year")
What kind of object did we create?
In [126]:
type(sl_grouped_year)
Out[126]:
What to do with that strange GroupBy object? We can first loop over it to get the labels and the sub-dataframes for each group:
In [127]:
for group_name, subdf in sl_grouped_year:
print(group_name)
print(subdf)
print("")
We could have done the same with less effort by grouping by the result of a custom function applied to the index. Let's reset the dataframe:
In [128]:
mean_sea_level = mean_sea_level.drop(["year"], axis=1).set_index("date")
So that we can do the groupby on the index:
In [129]:
sl_grouped_year = mean_sea_level.groupby(int)
Something else that can be done with such an object is to look at its groups attribute to see the labels mapped to the rows involved:
In [130]:
sl_grouped_year.groups
Out[130]:
How to aggregate the results of this grouping depends on what we want to see: do we want to see averaged over the years? That is so common that it has been implemented directly as a method on the GroupBy object.
In [131]:
sl_grouped_year.mean()
Out[131]:
In [132]:
# We can apply any other reduction function or even a dict of functions using aggregate:
sl_grouped_year.aggregate({"mean_global": np.std})
Out[132]:
Another possibility is to transform each group separately, rather than aggregate. For example, here we group over decades and subtract to each value, the average over that decade:
In [133]:
sl_grouped_decade = mean_sea_level.groupby(lambda x: int(x/10.))
sl_grouped_decade.groups.keys()
Out[133]:
In [134]:
sl_grouped_decade.transform(lambda subframe: (subframe - subframe.mean()/subframe.std()))
Out[134]:
Pivot table also allows to summarize the information, allowing to convert repeating columns into axes. For example, let's say that we would like to know how many sea level stations are in various european countries. And we would like to group the answers into 2 categories: the stations that have been updated recently (after 2000) and the others.
Let's first extract only entries located (roughly) in Europe.
In [135]:
european_filter = ((local_sea_level_stations["Lat"] > 30) &
(local_sea_level_stations["Lat"] < 70) &
(local_sea_level_stations["Lon"] > -10) &
(local_sea_level_stations["Lon"] < 40)
)
# Let's make a copy to work with a new, clean block of memory
# (if you are interested, try and remove the copy to see the consequences further down...
european_stations = local_sea_level_stations[european_filter].copy()
european_stations["Country"].unique()
Out[135]:
The columns of our future table should have 2 values, whether the station was updated recently or not. Let's build a column to store that information:
In [136]:
european_stations["Recently updated"] = european_stations["Date"] > pd.to_datetime("2000")
Finally, what value will be displayed inside the table. The values should be extracted from a column, pivot_table allowing an aggregation function to be applied when more than 1 value is found for a given case. Each station should count for 1, and we could aggregate multiple stations by summing these ones:
In [137]:
european_stations["Number of stations"] = np.ones(len(european_stations))
In [138]:
european_stations.sort_values(by="Country")
Out[138]:
In [139]:
station_counts = pd.pivot_table(european_stations, index="Country", columns="Recently updated",
values="Number of stations", aggfunc=np.sum)
# Let's remove from the table the countries for which no station was found:
station_counts.dropna(how="all")
Out[139]:
QUIZ: Why is there still some countries with no entries?
EXERCISE: How many recently updated stations? Not recently updated stations? Which country has the most stations? Which country has the most recently updated stations?
In [140]:
# Your code here
EXERCISE: How would we build the same dataframe with a groupby operation?
In [141]:
# Your code here
EXERCISE: Refer to exercises/pivot_table/pivot_tables.py
Both Series and dataframes have a corr method to compute the correlation coefficient between series:
In [142]:
# Let's see what how the various sea levels are correlated with each other:
mean_sea_level["northern_hem"].corr(mean_sea_level["southern_hem"])
Out[142]:
In [143]:
# If series are already grouped into a DataFrame, computing all correlation coeff is trivial:
mean_sea_level.corr()
Out[143]:
Note: by default, the method used is the Pearson correlation coefficient (https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient). Other methods are available (kendall, spearman using the method kwarg).
In [144]:
# Visualize the correlation matrix
plt.imshow(mean_sea_level.corr(), interpolation="nearest")
Out[144]:
In [145]:
plt.yticks?
In [146]:
# let's make it a little better to confirm that learning about global sea level cannot be done from just
# looking at stations in the northern hemisphere:
plt.imshow(mean_sea_level.corr(), interpolation="nearest")
plt.xticks(np.arange(3), mean_sea_level.corr().columns)
plt.yticks(np.arange(3), mean_sea_level.corr().index)
plt.colorbar()
Out[146]:
The recommeded way to build ordinaty least square regressions is by using statsmodels.
In [147]:
import statsmodels.formula.api as sm
In [148]:
sm_model = sm.ols(formula="mean_global ~ northern_hem + southern_hem", data=mean_sea_level).fit()
In [149]:
sm_model.params
Out[149]:
In [150]:
type(sm_model.params)
Out[150]:
In [151]:
sm_model.summary()
Out[151]:
In [152]:
plt.figure(figsize=LARGE_FIGSIZE)
mean_sea_level["mean_global"].plot()
sm_model.fittedvalues.plot(label="OLS prediction")
plt.legend(loc="upper left")
Out[152]:
Now, we would like to look for correlations between our monthly temperatures and the sea levels we have. For this to be possible, some data alignment must be done since the time scales are very different for the 2 datasets.
In [153]:
mean_sea_level["mean_global"].index
Out[153]:
In [154]:
giss_temp_series.index
Out[154]:
In [155]:
DAYS_PER_YEAR = {}
In [156]:
import calendar
# Let's first convert the floating point dates in the sea level to timestamps:
def floating_year_to_timestamp(float_date):
""" Convert a date as a floating point year number to a pandas timestamp object.
"""
year = int(float_date)
days_per_year = 366 if calendar.isleap(year) else 365
remainder = float_date - year
daynum = 1 + remainder * (days_per_year - 1)
daynum = int(round(daynum))
# Convert day number to month and day
day = daynum
month = 1
while month < 13:
month_days = calendar.monthrange(year, month)[1]
if day <= month_days:
return pd.Timestamp(str(year)+"/"+str(month)+"/"+str(day))
day -= month_days
month += 1
raise ValueError('{} does not have {} days'.format(year, daynum))
In [157]:
floating_year_to_timestamp(1996.0), floating_year_to_timestamp(1996.5), floating_year_to_timestamp(1996.9999)
Out[157]:
In [158]:
dt_index = pd.Series(mean_sea_level["mean_global"].index).apply(floating_year_to_timestamp)
dt_index
Out[158]:
In [159]:
mean_sea_level = mean_sea_level.reset_index(drop=True)
mean_sea_level.index = dt_index
mean_sea_level
Out[159]:
Now, how to align the 2 series? Is this one sampled regularly so that the month temperatures can be upscaled to that frequency?
What is the frequency of that new index?
In [160]:
dt_index.dtype
Out[160]:
In [161]:
# What is the frequency of the new index? The numpy way to compute differences between all values doesn't work:
dt_index[1:] - dt_index[:-1]
Out[161]:
IMPORTANT Note: The above failure is due to the fact that operations between series automatically align them based on their index.
In [162]:
# There is a method for shifting values up/down the index:
dt_index.shift()
Out[162]:
In [163]:
# So the distances can be computed with
dt_index - dt_index.shift()
Out[163]:
In [164]:
# Not constant reads apparently. Let's downscale the frequency of the sea levels
# to monthly, like the temperature reads we have:
monthly_mean_sea_level = mean_sea_level.resample("MS").mean().to_period()
monthly_mean_sea_level
Out[164]:
In [165]:
monthly_mean_sea_level["mean_global"].align(giss_temp_series)
Out[165]:
In [166]:
giss_temp_series.align?
In [167]:
# Now that the series are using the same type and frequency of indexes, to align them is trivial:
monthly_mean_sea_level["mean_global"].align(giss_temp_series, join='inner')
Out[167]:
In [168]:
aligned_sl, aligned_temp = monthly_mean_sea_level["mean_global"].align(giss_temp_series, join='inner')
aligned_df = pd.DataFrame({"mean_sea_level": aligned_sl, "mean_global_temp": aligned_temp})
The alignment can even be done on an entire dataframe:
In [169]:
monthly_mean_sea_level.align(giss_temp_series, axis=0, join='inner')
Out[169]:
In [170]:
aligned_sea_levels, aligned_temp = monthly_mean_sea_level.align(giss_temp_series, axis=0, join='inner')
aligned_monthly_data = aligned_sea_levels.copy()
aligned_monthly_data["global_temp"] = aligned_temp
aligned_monthly_data
Out[170]:
In [171]:
aligned_monthly_data.plot(figsize=LARGE_FIGSIZE)
Out[171]:
In [172]:
aligned_monthly_data.corr()
Out[172]:
In [173]:
model = sm.ols("southern_hem ~ global_temp", data=aligned_monthly_data)
params = model.fit()
params.rsquared
Out[173]:
What if we had done the analysis yearly instead of monthly to remove seasonal variations?
In [174]:
aligned_yearly_data = aligned_monthly_data.resample("A").mean()
aligned_yearly_data.plot()
Out[174]:
In [175]:
aligned_yearly_data.corr()
Out[175]:
In [176]:
model = sm.ols("southern_hem ~ global_temp", data=aligned_yearly_data).fit()
model.rsquared
Out[176]:
An auto-regresssive model fits existing data and build a (potentially predictive) model of the data fitted. We use the timeseries analysis (tsa) submodule of statsmodels to make out-of-sample predictions for the upcoming decades:
In [177]:
from statsmodels.tsa.api import AR
In [178]:
import statsmodels as sm
# Let's remove seasonal variations by resampling annually
data = giss_temp_series.resample("A").mean().to_timestamp()
ar_model = AR(data, freq='A')
ar_res = ar_model.fit(maxlag=60, disp=True)
In [179]:
plt.figure(figsize=LARGE_FIGSIZE)
pred = ar_res.predict(start='1950-1-1', end='2070')
data.plot(style='k', label="Historical Data")
pred.plot(style='r', label="Predicted Data")
plt.ylabel("Temperature variation (0.01 degC)")
plt.legend()
Out[179]:
EXERCISE: Make another auto-regression on the sea level of the Atlantic ocean to estimate how much New York is going to flood in the coming century.
data/sea_levels/sl_Atlantic_Ocean.txt. local_sea_level_stations dataset, and use it to download timeseries in NY (URL would be http://www.psmsl.org/data/obtaining/met.monthly.data/< ID >.metdata).
In [180]:
# Your code here