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: <p><font size="6"><b> Case study: air quality data of European monitoring stations (AirBase)</b></font></p>
<br>
AirBase (The European Air quality dataBase): hourly measurements of all air quality monitoring stations from Europe. 

© 2016, Joris Van den Bossche and Stijn Van Hoey  (&#106;&#111;&#114;&#105;&#115;&#118;&#97;&#110;&#100;&#101;&#110;&#98;&#111;&#115;&#115;&#99;&#104;&#101;&#64;&#103;&#109;&#97;&#105;&#108;&#46;&#99;&#111;&#109;, &#115;&#116;&#105;&#106;&#110;&#118;&#97;&#110;&#104;&#111;&#101;&#121;&#64;&#103;&#109;&#97;&#105;&#108;&#46;&#99;&#111;&#109;). Licensed under CC BY 4.0 Creative Commons


AirBase is the European air quality database maintained by the European Environment Agency (EEA). It contains air quality monitoring data and information submitted by participating countries throughout Europe. The air quality database consists of a multi-annual time series of air quality measurement data and statistics for a number of air pollutants.
End of explanation
"""


%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

pd.options.display.max_rows = 8

"""
Explanation: Some of the data files that are available from AirBase were included in the data folder: the hourly concentrations of nitrogen dioxide (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
End of explanation
"""


with open("data/BETR8010000800100hour.1-1-1990.31-12-2012") as f:
    print(f.readline())

"""
Explanation: Processing a single file
We will start with processing one of the downloaded files (BETR8010000800100hour.1-1-1990.31-12-2012). Looking at the data, you will see it does not look like a nice csv file:
End of explanation
"""


data = pd.read_csv("data/BETR8010000800100hour.1-1-1990.31-12-2012", sep='\t')#, header=None)

data.head()

"""
Explanation: So we will need to do some manual processing.
Just reading the tab-delimited data:
End of explanation
"""


# Column names: list consisting of 'date' and then intertwined the hour of the day and 'flag'
hours = ["{:02d}".format(i) for i in range(24)]
column_names = ['date'] + [item for pair in zip(hours, ['flag']*24) for item in pair]

# %load snippets/07 - Case study - air quality data7.py

# %load snippets/07 - Case study - air quality data8.py

"""
Explanation: The above data is clearly not ready to be used! Each row contains the 24 measurements for each hour of the day, and also contains a flag (0/1) indicating the quality of the data. Furthermore, there is no header row with column names.
<div class="alert alert-success">

<b>EXERCISE</b>: <br><br> Clean up this dataframe by using more options of `read_csv` (see its [docstring](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_csv.html))

 <ul>
  <li>specify the correct delimiter</li>
  <li>specify that the values of -999 and -9999 should be regarded as NaN</li>
  <li>specify are own column names (for how the column names are made up, see See http://stackoverflow.com/questions/6356041/python-intertwining-two-lists)
</ul>
</div>
End of explanation
"""


flag_columns = [col for col in data.columns if 'flag' in col]
# we can now use this list to drop these columns

# %load snippets/07 - Case study - air quality data10.py

data.head()

"""
Explanation: For the sake of this tutorial, we will disregard the 'flag' columns (indicating the quality of the data). 
<div class="alert alert-success">

<b>EXERCISE</b>:
<br><br>
Drop all 'flag' columns ('flag1', 'flag2', ...)
End of explanation
"""


# %load snippets/07 - Case study - air quality data12.py

# %load snippets/07 - Case study - air quality data13.py

# %load snippets/07 - Case study - air quality data14.py

# %load snippets/07 - Case study - air quality data15.py

# %load snippets/07 - Case study - air quality data16.py

data_stacked.head()

"""
Explanation: Now, we want to reshape it: our goal is to have the different hours as row indices, merged with the date into a datetime-index. Here we have a wide and long dataframe, and want to make this a long, narrow timeseries.
<div class="alert alert-info">

<b>REMEMBER</b>: 

 <ul>
  <li>Recap: reshaping your data with [`stack` and `unstack`](./pandas_07_reshaping_data.ipynb)</li>
</ul>


<img src="img/schema-stack.svg" width=70%>

</div>

<div class="alert alert-success">

<b>EXERCISE</b>:

<br><br>

Reshape the dataframe to a timeseries. 
The end result should look like:<br><br>


<div class='center'>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>BETR801</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>1990-01-02 09:00:00</th>
      <td>48.0</td>
    </tr>
    <tr>
      <th>1990-01-02 12:00:00</th>
      <td>48.0</td>
    </tr>
    <tr>
      <th>1990-01-02 13:00:00</th>
      <td>50.0</td>
    </tr>
    <tr>
      <th>1990-01-02 14:00:00</th>
      <td>55.0</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
    </tr>
    <tr>
      <th>2012-12-31 20:00:00</th>
      <td>16.5</td>
    </tr>
    <tr>
      <th>2012-12-31 21:00:00</th>
      <td>14.5</td>
    </tr>
    <tr>
      <th>2012-12-31 22:00:00</th>
      <td>16.5</td>
    </tr>
    <tr>
      <th>2012-12-31 23:00:00</th>
      <td>15.0</td>
    </tr>
  </tbody>
</table>
<p style="text-align:center">170794 rows × 1 columns</p>
</div>

 <ul>
  <li>Reshape the dataframe so that each row consists of one observation for one date + hour combination</li>
  <li>When you have the date and hour values as two columns, combine these columns into a datetime (tip: string columns can be summed to concatenate the strings) and remove the original columns</li>
  <li>Set the new datetime values as the index, and remove the original columns with date and hour values</li>

</ul>


**NOTE**: This is an advanced exercise. Do not spend too much time on it and don't hesitate to look at the solutions. 

</div>
End of explanation
"""


data_stacked.index

data_stacked.plot()

"""
Explanation: Our final data is now a time series. In pandas, this means that the index is a DatetimeIndex:
End of explanation
"""


def read_airbase_file(filename, station):
    """
    Read hourly AirBase data files.
    
    Parameters
    ----------
    filename : string
        Path to the data file.
    station : string
        Name of the station.
       
    Returns
    -------
    DataFrame
        Processed dataframe.
    """
    
    ...
    
    return ...

# %load snippets/07 - Case study - air quality data21.py

"""
Explanation: Processing a collection of files
We now have seen the code steps to process one of the files. We have however multiple files for the different stations with the same structure. Therefore, to not have to repeat the actual code, let's make a function from the steps we have seen above.
<div class="alert alert-success">

<b>EXERCISE</b>:

 <ul>
  <li>Write a function `read_airbase_file(filename, station)`, using the above steps the read in and process the data, and that returns a processed timeseries.</li>
</ul>
</div>
End of explanation
"""


filename = "data/BETR8010000800100hour.1-1-1990.31-12-2012"
station = filename.split("/")[-1][:7]

station

test = read_airbase_file(filename, station)
test.head()

"""
Explanation: Test the function on the data file from above:
End of explanation
"""


import glob

# %load snippets/07 - Case study - air quality data33.py

"""
Explanation: We now want to use this function to read in all the different data files from AirBase, and combine them into one Dataframe. 
<div class="alert alert-success">

<b>EXERCISE</b>:

 <ul>
  <li>Use the `glob.glob` function to list all 4 AirBase data files that are included in the 'data' directory, and call the result `data_files`.</li>
</ul>
</div>
End of explanation
"""


# %load snippets/07 - Case study - air quality data34.py

# %load snippets/07 - Case study - air quality data35.py

combined_data.head()

"""
Explanation: <div class="alert alert-success">

<b>EXERCISE</b>:

 <ul>
  <li>Loop over the data files, read and process the file using our defined function, and append the dataframe to a list.</li>
  <li>Combine the the different DataFrames in the list into a single DataFrame where the different columns are the different stations. Call the result `combined_data`.</li>

</ul>
</div>
End of explanation
"""


combined_data.to_csv("airbase_data.csv")

"""
Explanation: Finally, we don't want to have to repeat this each time we use the data. Therefore, let's save the processed data to a csv file.
End of explanation
"""


alldata = pd.read_csv('airbase_data.csv', index_col=0, parse_dates=True)

"""
Explanation: Working with time series data
We processed the individual data files above, and saved it to a csv file airbase_data.csv. Let's import the file here (if you didn't finish the above exercises, a version of the dataset is also available in data/airbase_data.csv):
End of explanation
"""


data = alldata['1999':].copy()

"""
Explanation: We only use the data from 1999 onwards:
End of explanation
"""


data.head() # tail()

data.info()

data.describe(percentiles=[0.1, 0.5, 0.9])

"""
Explanation: Som first exploration with the typical functions:
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
"""


# %load snippets/07 - Case study - air quality data50.py

# %load snippets/07 - Case study - air quality data51.py

"""
Explanation: Exercises
<div class="alert alert-warning">

<b>REMINDER</b>: <br><br>

Take a look at the [Timeseries notebook](05 - Time series data.ipynb) when you require more info about...

 <ul>
  <li>`resample`</li>
  <li>string indexing of DateTimeIndex</li>
</ul><br><br>

</div>

<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
"""


# %load snippets/07 - Case study - air quality data52.py

# %load snippets/07 - Case study - air quality data53.py

"""
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
"""


# %load snippets/07 - Case study - air quality data54.py

"""
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
"""


# %load snippets/07 - Case study - air quality data55.py

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: The evolution of the yearly averages with, and the overall mean of all stations (indicate the overall mean with a thicker black line)?
</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
"""


# %load snippets/07 - Case study - air quality data57.py

"""
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
"""


# %load snippets/07 - Case study - air quality data58.py

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


# %load snippets/07 - Case study - air quality data59.py

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

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


df2011 = data['2011'].dropna()

# %load snippets/07 - Case study - air quality data64.py

# %load snippets/07 - Case study - air quality data65.py

"""
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
"""


# %load snippets/07 - Case study - air quality data66.py

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


# %load snippets/07 - Case study - air quality data67.py

# %load snippets/07 - Case study - air quality data68.py

# %load snippets/07 - Case study - air quality data69.py

# %load snippets/07 - Case study - air quality data70.py

"""
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 for each year/station?
</div>
End of explanation
"""


# %load snippets/07 - Case study - air quality data72.py

# %load snippets/07 - Case study - air quality data73.py

# %load snippets/07 - Case study - air quality data74.py

"""
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
"""


data.index.weekday?

# %load snippets/07 - Case study - air quality data76.py

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


# %load snippets/07 - Case study - air quality data77.py

# %load snippets/07 - Case study - air quality data78.py

# %load snippets/07 - Case study - air quality data79.py

# %load snippets/07 - Case study - air quality data80.py

# %load snippets/07 - Case study - air quality data81.py

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


# %load snippets/07 - Case study - air quality data82.py

# %load snippets/07 - Case study - air quality data83.py

# %load snippets/07 - Case study - air quality data84.py

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: Visualize the typical week profile for the different stations as boxplots (where the values in one boxplot are the daily means for the different weeks for a certain weekday).
</div>

Tip: the boxplot method of a DataFrame expects the data for the different boxes in different columns). For this, you can either use pivot_table as a combination of groupby and unstack
End of explanation
"""


# %load snippets/07 - Case study - air quality data85.py

"""
Explanation: An alternative method using groupby and unstack:
End of explanation
"""


# %load snippets/07 - Case study - air quality data86.py

"""
Explanation: <div class="alert alert-success">
    <b>QUESTION</b>: The maximum daily 8 hour mean should be below 100 µg/m³. What are the number of exceedances of this limit for each year/station?
</div>

Tip: have a look at the rolling method to perform moving window operations.
Note: this is not an actual limit for NO2, but a nice exercise to introduce the rolling method. Other pollutans, such as 03 have actually such kind of limit values.
End of explanation
"""


# %load snippets/07 - Case study - air quality data87.py

# %load snippets/07 - Case study - air quality data88.py

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


