%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
try:
    import seaborn
except:
    pass

pd.options.display.max_rows = 8

"""
Explanation: Working with time series data
Some imports:
End of explanation
"""


from IPython.display import HTML
HTML('<iframe src=http://www.eea.europa.eu/data-and-maps/data/airbase-the-european-air-quality-database-8#tab-data-by-country width=900 height=350></iframe>')

"""
Explanation: Case study: air quality data of European monitoring stations (AirBase)
AirBase (The European Air quality dataBase): hourly measurements of all air quality monitoring stations from Europe.
End of explanation
"""


!head -5 data/airbase_data.csv

"""
Explanation: I downloaded and preprocessed some of the data (python-airbase): data/airbase_data.csv. This file includes the hourly concentrations of NO2 for 4 different measurement stations:

FR04037 (PARIS 13eme): urban background site at Square de Choisy
FR04012 (Paris, Place Victor Basch): urban traffic site at Rue d'Alesia
BETR802: urban traffic site in Antwerp, Belgium
BETN029: rural background site in Houtem, Belgium

See http://www.eea.europa.eu/themes/air/interactive/no2
Importing the data
Import the csv file:
End of explanation
"""


data = pd.read_csv('data/airbase_data.csv', index_col=0, parse_dates=True, na_values=[-9999])

"""
Explanation: As you can see, the missing values are indicated by -9999. This can be recognized by read_csv by passing the na_values keyword:
End of explanation
"""


data.head(3)

data.tail()

"""
Explanation: Exploring the data
Some useful methods:
head and tail
End of explanation
"""


data.info()

"""
Explanation: info()
End of explanation
"""


data.describe()

"""
Explanation: Getting some basic summary statistics about the data with describe:
End of explanation
"""


data.plot(kind='box', ylim=[0,250])

data['BETR801'].plot(kind='hist', bins=50)

data.plot(figsize=(12,6))

"""
Explanation: Quickly visualizing the data
End of explanation
"""


data[-500:].plot(figsize=(12,6))

"""
Explanation: This does not say too much ..
We can select part of the data (eg the latest 500 data points):
End of explanation
"""


data.index

"""
Explanation: Or we can use some more advanced time series features -> next section!
Working with time series data
When we ensure the DataFrame has a DatetimeIndex, time-series related functionality becomes available:
End of explanation
"""


data["2010-01-01 09:00": "2010-01-01 12:00"]

"""
Explanation: Indexing a time series works with strings:
End of explanation
"""


data['2012']

"""
Explanation: A nice feature is "partial string" indexing, where we can do implicit slicing by providing a partial datetime string.
E.g. all data of 2012:
End of explanation
"""


data['2012-01':'2012-03']

"""
Explanation: Normally you would expect this to access a column named '2012', but as for a DatetimeIndex, pandas also tries to interprete it as a datetime slice.
Or all data of January up to March 2012:
End of explanation
"""


data.index.hour

data.index.year

"""
Explanation: Time and date components can be accessed from the index:
End of explanation
"""


data = data['1999':]

"""
Explanation: <div class="alert alert-success">
    <b>EXERCISE</b>: select all data starting from 1999
</div>
End of explanation
"""


data[data.index.month == 1]

"""
Explanation: <div class="alert alert-success">
    <b>EXERCISE</b>: select all data in January for all different years
</div>
End of explanation
"""


data['months'] = data.index.month
data[data['months'].isin([1, 2, 3])]

"""
Explanation: <div class="alert alert-success">
    <b>EXERCISE</b>: select all data in January, February and March for all different years
</div>
End of explanation
"""


data[(data.index.hour >= 8) & (data.index.hour < 20)]

data.between_time('08:00', '20:00')

"""
Explanation: <div class="alert alert-success">
    <b>EXERCISE</b>: select all 'daytime' data (between 8h and 20h) for all days
</div>
End of explanation
"""


data.resample('D').head()

"""
Explanation: The power of pandas: resample
A very powerfull method is resample: converting the frequency of the time series (e.g. from hourly to daily data).
The time series has a frequency of 1 hour. I want to change this to daily:
End of explanation
"""


data.resample('D', how='max').head()

"""
Explanation: By default, resample takes the mean as aggregation function, but other methods can also be specified:
End of explanation
"""


data.resample('M').plot() # 'A'

# data['2012'].resample('D').plot()

"""
Explanation: The string to specify the new time frequency: http://pandas.pydata.org/pandas-docs/dev/timeseries.html#offset-aliases
These strings can also be combined with numbers, eg '10D'.
Further exploring the data:
End of explanation
"""


data.loc['2009':, 'FR04037'].resample('M', how='mean').plot()
data.loc['2009':, 'FR04037'].resample('M', how='median').plot()

data.loc['2009':, 'FR04037'].resample('M', how=['mean', 'median']).plot()

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: plot the monthly mean and median concentration of the 'FR04037' station for the years 2009-2012
</div>
End of explanation
"""


daily = data['FR04037'].resample('D')

daily.resample('M', how=['min', 'max']).plot()

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: plot the monthly mininum and maximum daily concentration of the 'BETR801' station
</div>
End of explanation
"""


data['2012'].mean().plot(kind='bar')

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: make a bar plot of the mean of the stations in year of 2012
</div>
End of explanation
"""


data.resample('A').plot()
data.mean(axis=1).resample('A').plot(color='k', linestyle='--', linewidth=4)

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: The evolution of the yearly averages with, and the overall mean of all stations?
</div>
End of explanation
"""


data.groupby(data.index.year).mean().plot()

"""
Explanation: Combination with groupby
resample can actually be seen as a specific kind of groupby. E.g. taking annual means with data.resample('A', 'mean') is equivalent to data.groupby(data.index.year).mean() (only the result of resample still has a DatetimeIndex).
End of explanation
"""


data['month'] = data.index.month

"""
Explanation: But, groupby is more flexible and can also do resamples that do not result in a new continuous time series, e.g. by grouping by the hour of the day to get the diurnal cycle.
<div class="alert alert-success">
    <b>QUESTION</b>: how does the *typical monthly profile* look like for the different stations?
</div>

1. add a column to the dataframe that indicates the month (integer value of 1 to 12):
End of explanation
"""


data.groupby('month').mean()

"""
Explanation: 2. Now, we can calculate the mean of each month over the different years:
End of explanation
"""


data.groupby('month').mean().plot()

"""
Explanation: 3. plot the typical monthly profile of the different stations:
End of explanation
"""


df2011 = data['2011']
df2011.groupby(df2011.index.week)[['BETN029', 'BETR801']].quantile(0.95).plot()

data = data.drop('month', axis=1)

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: plot the weekly 95% percentiles of the concentration in 'BETR801' and 'BETN029' for 2011
</div>
End of explanation
"""


data.groupby(data.index.hour).mean().plot()

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: The typical diurnal profile for the different stations?
</div>
End of explanation
"""


data.index.weekday?

data['weekday'] = data.index.weekday

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: What is the difference in the typical diurnal profile between week and weekend days?
</div>
End of explanation
"""


data['weekend'] = data['weekday'].isin([5, 6])

data_weekend = data.groupby(['weekend', data.index.hour]).mean()
data_weekend.head()

data_weekend_FR04012 = data_weekend['FR04012'].unstack(level=0)
data_weekend_FR04012.head()

data_weekend_FR04012.plot()

data = data.drop(['weekday', 'weekend'], axis=1)

"""
Explanation: Add a column indicating week/weekend
End of explanation
"""


exceedances = data > 200

# group by year and count exceedances (sum of boolean)
exceedances = exceedances.groupby(exceedances.index.year).sum()

exceedances

ax = exceedances.loc[2005:].plot(kind='bar')
ax.axhline(18, color='k', linestyle='--')

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: What are the number of exceedances of hourly values above the European limit 200 µg/m3 ?
</div>
End of explanation
"""


yearly = data['2000':].resample('A')

(yearly > 40).sum()

yearly.plot()
plt.axhline(40, linestyle='--', color='k')

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: And are there exceedances of the yearly limit value of 40 µg/m3 since 200 ?
</div>
End of explanation
"""


# add a weekday and week column
data['weekday'] = data.index.weekday
data['week'] = data.index.week
data.head()

# pivot table so that the weekdays are the different columns
data_pivoted = data['2012'].pivot_table(columns='weekday', index='week', values='FR04037')
data_pivoted.head()

box = data_pivoted.boxplot()

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: Visualize the typical week profile for the different stations as boxplots.
</div>

Tip: the boxplot method of a DataFrame expects the data for the different boxes in different columns)
End of explanation
"""


data[['BETR801', 'BETN029', 'FR04037', 'FR04012']].corr()

data[['BETR801', 'BETN029', 'FR04037', 'FR04012']].resample('D').corr()

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: Calculate the correlation between the different stations
</div>
End of explanation
"""


