forecasting - freeCodeCamp.org

How I Built My Own Forecasting Tool Using a Weather API

freeCodeCamp — Mon, 11 May 2020 17:50:21 +0000

By Anton Lawrence

Abrupt weather and climate change are things that everybody is dealing with. In fact, the vast majority of the global population relies on accurate, real-time weather data and forecasts to make informed decisions.

This has increased the importance of reliable Android and iOS weather apps. In this article, we will show you how to create a simple forecasting tool using NodeJS and a weather API.

But before that, let’s go over the importance of weather apps.

Why Do We Need Weather Apps?

A feature-rich weather forecasting app can provide great value to various industries. Some noteworthy benefits of weather applications include:

Providing immediate access to local weather conditions and upcoming forecast, thus saving you time.
Providing real-time notifications about prevailing and expected weather conditions.
Helping governments and local administrations prepare for natural disasters and save lives.
Helping farmers take preventive measures.
Facilitating the global travel and tourism industry.
Providing clear weather forecasts which are crucial in the aviation and logistics industry.

What You Need to Build a Weather App

Here are some of the things you will need to build a weather app successfully:

Familiarity working with JavaScript (Node.js)
A text editor such as Notepad or an IDE. My favorite is Visual Studio.
Access to a reliable weather API such as ClimaCell
Access to a map service
Knowledge of HTML, CSS, and Bootstrap

Once you have these ready, you’re good to go.

Overview of the ClimaCell Weather API

ClimaCell is a popular weather provider that offers hyper-accurate historical weather data as well as forecasts through an easy to consume API.

The Building Process

In this section, I will show you how I created a forecasting app where a user enters their city or any other location by name and fetches weather data from the ClimaCell API. The API responds to the request by returning data, which is then displayed to the user.

Install NodeJS and Create a New Project

For this project, we will use Node.js — one of the most popular run-time environments for JavaScript. Node.js helps developers create quick web applications. It has a wide range of libraries and modules for creating advanced web applications.

If you do not have Node.js on your device, you can install it from the official website.
Once installed, we use this command to initialize npm - the default packet manager used by Node.js.

$ npm init

This creates our project, so you will be prompted to enter a few details such as package name, description, Git repository, and more.

Next, we install the modules required to run our project. To generate a Node.js app skeleton, we use express - a framework for building Node.js web applications.

$ npm install express

Installing the express framework helps you run the server, handle client requests, and connect the right HTML template with a response.

Next, we will also install unirest - a simple yet powerful solution that allows you to request a library.

It will help us make requests to the ClimaCell API and handle the responses.

Use this command:

npm install unirest

At this point, we’ve installed the necessary modules, and the project is ready.
Next, we generate a weather app using the express generator tool. On the command line, type this:

express --view=pug weather-app-nodejs

You should now have a view like this on the command line:

Get the ClimaCell Weather API

To get access to the ClimaCell API, you will need to sign up for an account on their page.

Once you create an account, sign in to their Microweather API dashboard which looks like this:

On the dashboard, click on references to check the API endpoints. As you can see, the ClimaCell API has a number of endpoints including the short term forecast, hourly forecast, real time data, and more.

Worth a mention is that each endpoint has its own code snippet. For example, here is the Node.js code snippet for getting real time weather data.

[Raw](https://carbon.now.sh/embed?bg=rgba(0%2C0%2C0%2C1)&t=blackboard&wt=none&l=javascript&ds=true&dsyoff=20px&dsblur=68px&wc=true&wa=true&pv=56px&ph=56px&ln=false&fl=1&fm=Hack&fs=14px&lh=133%25&si=false&es=1x&wm=false&code=router.post(%27%252Fweather%27%252C%2520function(req%252C%2520res%252C%2520next)%257B%250A%250A%2520%2520let%2520city%2520%253D%2520req.body.city%253B%250A%250A%2520%2520url%2520%253D%2520url%252Bcity%252B%2522%2526%2522%252BappId%253B%250A%250A%2520request(url%252C%2520function%2520(error%252C%2520response%252C%2520body)%2520%257B%250A%250A%2520%2520%2520%2520%2520%2520body%2520%253D%2520JSON.parse(body)%253B%250A%250A%2520%2520%2520%2520%2520%2520if(error%2520%2526%2526%2520response.statusCode%2520!%253D%2520200)%257B%250A%250A%2520%2520%2520%2520%2520%2520%2520%2520throw%2520error%253B%250A%250A%2520%2520%2520%2520%2520%2520%257D%250A%250A%2520%2520%2520%2520let%2520country%2520%253D%2520(body.sys.country)%2520%253F%2520body.sys.country%2520%253A%2520%27%27%2520%253B%250A%250A%2520%2520%2520%2520let%2520var%2520request%2520%253D%2520require(%2522request%2522)%253B%250A%250Avar%2520options%2520%253D%2520%257B%250A%2520%2520method%253A%2520%27GET%27%252C%250A%2520%2520url%253A%2520%27https%253A%252F%252Fapi.climacell.co%252Fv3%252Fweather%252Frealtime%27%252C%250A%2520%2520qs%253A%2520%257Bapikey%253A%2520%279efRx8KrGKa6ME0ORWIJz7TaNjuRQAvb%27%257D%250A%257D%253B%250A%250Arequest(options%252C%2520function%2520(error%252C%2520response%252C%2520body)%2520%257B%250A%2520%2520if%2520(error)%2520throw%2520new%2520Error(error)%253B%250A%250A%2520%2520console.log(body)%253B%250A%257D)%253B%2520%253D%2520%2522For%2520city%2520%2522%252Bcity%252B%27%252C%2520country%2520%27%252Bcountry%253B%250A%250A%2520%2520%2520%2520res.render(%27index%27%252C%2520%257Bbody%2520%253A%2520body%252C%2520forecast%253A%2520forecast%257D)%253B%250A%250A%2520%2520%2520%257D)%253B%250A%250A%257D)%253B" rel="noreferrer nofollow noopener)

Modifying the application

To call the ClimaCell API, we need to first edit some files. Here, you can use Notepad or open the project directory in your IDE for easier editing. It should appear as shown:

We start modifying our files by adding bootstrap to layout.pug. Open the views directory and insert this snippet to the file.

Next, we create a form by adding the snippet below to index.pug file.

Notice how we use the HTTP post method to send data to the server. The code above also sets the action parameter to weather route and adds the text input as “city.”
An input button to fetch the weather is also added.

We now create an HTML table just below the form to display fetched weather records.

[Raw](https://carbon.now.sh/embed?bg=rgba(0%2C0%2C0%2C1)&t=blackboard&wt=none&l=htmlmixed&ds=true&dsyoff=20px&dsblur=68px&wc=true&wa=true&pv=56px&ph=56px&ln=false&fl=1&fm=Hack&fs=14px&lh=133%25&si=false&es=1x&wm=false&code=.row%250A%2520%2520%2520%2520%2520%2520%2520.col-md-12%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520p%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520br%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520br%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520strong%2520Weather%2520Forecast%2520%2523%257Bforecast%257D%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520table.table%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520thead%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520tr%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520th%2520Longitude%2520%252F%2520Latitude%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520th%2520Pressure%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520th%2520Temprature%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520th%2520Humidty%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520tbody%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520if%2520body%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520tr%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520td%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2523%257Bbody.coord.lon%257D%2520%252F%2520%2523%257Bbody.coord.lon%257D%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520td%2520%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2523%257Bbody.main.pressure%257D%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520td%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2523%257Bbody.main.temp%257D%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520td%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2523%257Bbody.main.humidity%257D%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520else%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520tr%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520td(colspan%253D%25226%2522)%250A%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%2520%257C%2520Enter%2520city%2520name%2520and%2520click%2520Fetch%2520Weather%2520button" rel="noreferrer nofollow noopener)

Inserting the code snippet above creates a table that looks like this:

Calling the ClimaCell API

To send requests to the ClimaCell API, we must install the request module.

npm i request --save

Next, we add the ClimaCell API credentials in the index.js file. Open the file in your routes directory and add the API key you obtained on the ClimaCell dashboard:

Here is the code to add API credentials:

After adding the API credentials, we update the index route. This is done by replacing the code section in ‘/’ route in index.js file.

We finish by creating the weather route in index.js.

[Raw](https://carbon.now.sh/embed?bg=rgba(0%2C0%2C0%2C1)&t=blackboard&wt=none&l=javascript&ds=true&dsyoff=20px&dsblur=68px&wc=true&wa=true&pv=56px&ph=56px&ln=false&fl=1&fm=Hack&fs=14px&lh=133%25&si=false&es=1x&wm=false&code=router.post(%27%252Fweather%27%252C%2520function(req%252C%2520res%252C%2520next)%257B%250A%250A%2520%2520let%2520city%2520%253D%2520req.body.city%253B%250A%250A%2520%2520url%2520%253D%2520url%252Bcity%252B%2522%2526%2522%252BappId%253B%250A%250A%2520request(url%252C%2520function%2520(error%252C%2520response%252C%2520body)%2520%257B%250A%250A%2520%2520%2520%2520%2520%2520body%2520%253D%2520JSON.parse(body)%253B%250A%250A%2520%2520%2520%2520%2520%2520if(error%2520%2526%2526%2520response.statusCode%2520!%253D%2520200)%257B%250A%250A%2520%2520%2520%2520%2520%2520%2520%2520throw%2520error%253B%250A%250A%2520%2520%2520%2520%2520%2520%257D%250A%250A%2520%2520%2520%2520let%2520country%2520%253D%2520(body.sys.country)%2520%253F%2520body.sys.country%2520%253A%2520%27%27%2520%253B%250A%250A%2520%2520%2520%2520let%2520forecast%2520%253D%2520%2522For%2520city%2520%2522%252Bcity%252B%27%252C%2520country%2520%27%252Bcountry%253B%250A%250A%2520%2520%2520%2520res.render(%27index%27%252C%2520%257Bbody%2520%253A%2520body%252C%2520forecast%253A%2520forecast%257D)%253B%250A%250A%2520%2520%2520%257D)%253B%250A%250A%257D)%253B" rel="noreferrer nofollow noopener)

This code snippet enables data in the input form to be posted to the index route. Once the user enters a city name, it is assigned to the city variable using the request object.

The URL is then appended with the city name and ID and the request sent to ClimaCell API.

The ClimaCell API server response is returned as a JSON file, which is then parsed and fed to the output template.

For instance, if the user was looking for Boston weather forecasts, the app will return this:

Note - The temperature in this example is shown in Kelvin and is equal to 50°F or 10°C.

Add Maps to Visualize Your Data

You can integrate interactive maps into your forecasting application to enhance user experience. This can be achieved by using a third-party map service provider for web applications.

Mapbox is one such tool that helps developers create awesome weather maps for their applications. It integrates seamlessly with any weather app.

To use the Mapbox, sign up on their website and check out their API. There are integrations for Android, iOS, Web, and Unity. In this case, we choose Web integration for our tool.

We can either install the Mapbox CDN or use the module bundler. Let's use the module bundler.

The first step would be installing the package

npm install Mapbox-gl –save

Next, we add the GL JS CSS file in the HTML file by including this snippet in the

We can now add the map to our application. To do this, use the code snippet below.

_[Raw](https://carbon.now.sh/embed?bg=rgba(0%2C0%2C0%2C1)&t=blackboard&wt=none&l=javascript&ds=true&dsyoff=20px&dsblur=68px&wc=true&wa=true&pv=56px&ph=56px&ln=false&fl=1&fm=Hack&fs=14px&lh=133%25&si=false&es=1x&wm=false&code=var%2520mapboxgl%2520%253D%2520require(%27mapbox-gl%252Fdist%252Fmapbox-gl.js%27)%253B%250A%2520%250Amapboxgl.accessToken%2520%253D%2520%27pk.eyJ1IjoiZGlja3Nvbi1tIiwiYSI6ImNrOXphd3MzZDBlMXYzbHFwM2kwbmlvbmkifQ.VE3RRbb8l78w9kxfmh_9ew%27%253B%250Avar%2520map%2520%253D%2520new%2520mapboxgl.Map(%257B%250Acontainer%253A%2520%27CONTAINER_ELEMENTID%27%252C%250Astyle%253A%2520%27mapbox%253A%252F%252Fstyles%252Fmapbox%252Fstreets-v11%27%250A%257D)%253B%250A" rel="noreferrer nofollow noopener)

You can choose where to place the map by replacing the

“CONTAINER_ELEMENT_ID’ .

Here is a sample map generated using Mapbox:

What’s next?

At this point, much of the work is done, and your app can get weather forecasts for any city using the ClimaCell API.

However, you can consider adding more interactive features to your application or extending its functionality.

Here are some things you might want to do:

Add a search function.
Improve the look of your user interface.
Query the application by ID or name.
Display a list of target cities and their respective IDs.
Add parameters to display additional weather data .
Integrate real-time notifications and warning signals.

As you can see, the basic app building process is pretty simple and straightforward. By following the above process to leverage the power of a weather API, even beginner-level developers can get their weather application up and running in a matter of minutes.

How to Forecast Air Pollution with Recurrent Neural Networks

freeCodeCamp — Thu, 02 Jan 2020 21:52:42 +0000

By Bert Carremans

After the citizen science project of Curieuze Neuzen, I wanted to learn more about air pollution to see if I could make a data science project out of it. On the website of the European Environment Agency, you can find a huge amount of data and information about air pollution.

In this notebook, we will focus on the air quality in Belgium, and more specifically on the pollution by sulphur dioxide (SO2). The data can be downloaded via https://www.eea.europa.eu/data-and-maps/data/aqereporting-2/be.

The zip file contains separate files for different air pollutants and aggregation levels. The first digit represents the pollutant ID as described in the vocabulary. The file used in this notebook is BE_1_2013–2015_aggregated_timeseries.csv. This is the SO2 pollution in Belgium, but you can also find similar data for other European countries.

Descriptions of the fields in the CSV files are available on the data download page. More background information on air pollutants can be found on Wikipedia.

Project Set-up

# Importing packages
from pathlib import Path
import pandas as pd
import numpy as np
import pandas_profiling
%matplotlib inline
import matplotlib.pyplot as plt
import warnings
warnings.simplefilter(action = 'ignore', category = FutureWarning)
from sklearn.preprocessing import MinMaxScaler

from keras.preprocessing.sequence import TimeseriesGenerator
from keras.models import Sequential
from keras.layers import Dense, LSTM, SimpleRNN
from keras.optimizers import RMSprop
from keras.callbacks import ModelCheckpoint, EarlyStopping
from keras.models import model_from_json

# Setting the project directory
project_dir = Path('/Users/bertcarremans/Data Science/Projecten/air_pollution_forecasting')

Loading the data

date_vars = ['DatetimeBegin','DatetimeEnd']

agg_ts = pd.read_csv(project_dir / 'data/raw/BE_1_2013-2015_aggregated_timeseries.csv', sep='\t', parse_dates=date_vars, date_parser=pd.to_datetime)
meta = pd.read_csv(project_dir / 'data/raw/BE_2013-2015_metadata.csv', sep='\t')

print('aggregated timeseries shape:{}'.format(agg_ts.shape))
print('metadata shape:{}'.format(meta.shape))

Data Exploration

Let’s use pandas_profiling to inspect the data.

pandas_profiling.ProfileReport(agg_ts)

I won’t show the output of pandas_profiling in this story in order not to clutter it with charts. But you can find it in my GitHub repo.

The pandas_profiling report shows us the following:

There are 6 constant variables. We can remove these from the data set.
No missing values exist, so probably we will not need to apply imputation.
AirPollutionLevel has some zeroes, but this could be perfectly normal. On the other hand, these variables have some extreme values, which might be incorrect recordings of air pollution.
There are 53 AirQualityStations, which are probably the same as the SamplingPoints. AirQualityStationEoICode is simply a shorter code for the AirQualityStation, so that variable can also be removed.
There are 3 values for AirQualityNetwork (Brussels, Flanders and Wallonia). Most measurements come from Flanders.
DataAggregationProcess: most rows contain data aggregated as the 24-hour mean of one day of measurements (P1D). More information on the other values can be found here. In this project, we will only consider P1D values.
DataCapture: Proportion of valid measurement time relative to the total measured time (time coverage) in the averaging period, expressed as a percentage. Almost all rows have about 100% of valid measurement time. Some rows have a DataCapture that is slightly lower than 100%.
DataCoverage: Proportion of valid measurement included in the aggregation process within the averaging period, expressed as a percentage. In this data set, we have a minimum of 75%. According to the definition of this variable values below 75% should not be included for air quality assessments, which explains why these rows are not present in the data set.
TimeCoverage: highly correlated to DataCoverage and will be removed from the data.
UnitOfAirPollutionLevel: 423 rows have a unit of count. To have a consistent target variable we will remove the records with this type of unit.
DateTimeBegin and DateTimeEnd: the histogram does not provide enough detail here. This needs to be analyzed further.

DateTimeBegin and DateTimeEnd

The histogram in the pandas_profiling combined multiple days per bin. Let’s look at a daily level how these variables behave.

Multiple aggregation levels per date

DatetimeBegin: a large number of records on the 1st of January of 2013, 2014, 2015 and 1st of October of 2013 and 2014.
DatetimeEnd: a large number of records on the 1st of January of 2014, 2015, 2016 and 1st of April of 2014 and 2015.

plt.figure(figsize=(20,6))
plt.plot(agg_ts.groupby('DatetimeBegin').count(), 'o', color='skyblue')
plt.title('Nb of measurements per DatetimeBegin')
plt.ylabel('number of measurements')
plt.xlabel('DatetimeBegin')
plt.show()

Number of rows per date

The outliers in the number of records are related to the multiple aggregation levels (DataAggregationProcess). The values in DataAggregationProcess on these dates reflect the time period between DatetimeBegin and DatetimeEnd. For example, the 1st of January 2013 is the start date of a one-year measurement period until the 1st of January 2014.

As we are only interested in the daily aggregation level, filtering out the other aggregation levels will solve this issue. We can also remove DatetimeEnd for that reason.

Missing timesteps on daily aggregation level

As we can see below, not all SamplingPoints have data for all DatetimeBegin in the three-year period. These are most likely days where the DataCoverage variable was below 75%. So on these days, we do not have sufficient valid measurements. Later in this notebook, we will use measurements on prior days to predict the pollution on the current day.

To have similarly sized timesteps, we will need to insert rows for the missing DatetimeBegin per SamplingPoint. We will insert the measurement data of the next day with valid data.

Secondly, we will remove the SamplingPoints with too many missing timesteps. Here we will take an arbitrary number of 1.000 timesteps as the minimum number of required timesteps.

ser_avail_days = agg_ts.groupby('SamplingPoint').nunique()['DatetimeBegin']
plt.figure(figsize=(8,4))
plt.hist(ser_avail_days.sort_values(ascending=False))
plt.ylabel('Nb SamplingPoints')
plt.xlabel('Nb of Unique DatetimeBegin')
plt.title('Distribution of Samplingpoints by the Nb of available measurement days')
plt.show()

Distribution of SamplingPoints by the number of available measurement days

Data Preparation

Data Cleaning

Based on the data exploration, we will do the following to clean the data:

Keeping only records with DataAggregationProcess of P1D
Removing records with UnitOfAirPollutionLevel of count
Removing unary variables and other redundant variables
Removing SamplingPoints which have less than 1000 measurement days

df = agg_ts.loc[agg_ts.DataAggregationProcess=='P1D', :] 
df = df.loc[df.UnitOfAirPollutionLevel!='count', :]
df = df.loc[df.SamplingPoint.isin(ser_avail_days[ser_avail_days.values >= 1000].index), :]
vars_to_drop = ['AirPollutant','AirPollutantCode','Countrycode','Namespace','TimeCoverage','Validity','Verification','AirQualityStation',
               'AirQualityStationEoICode','DataAggregationProcess','UnitOfAirPollutionLevel', 'DatetimeEnd', 'AirQualityNetwork',
               'DataCapture', 'DataCoverage']
df.drop(columns=vars_to_drop, axis='columns', inplace=True)

Inserting rows for the missing timesteps

For each SamplingPoint, we will first insert (empty) rows for which we do not have a DatetimeBegin. This can be done by creating a complete multi-index with all SamplingPoints and over the range between the minimum and maximum DatetimeBegin. Then, reindex will insert the missing rows but with NaN for the columns.

Secondly, we use bfill and specify to impute the missing values with the values of the next row with valid data. The bfill method is applied to a groupby object to limit the backfilling within the rows of each SamplingPoint. That way we do not use the values of another SamplingPoint to fill in the missing values.

A samplepoint to test whether this operation worked correctly is _SPO-BETR223_00001100 for the date 2013–01–29.

dates = list(pd.period_range(min(df.DatetimeBegin), max(df.DatetimeBegin), freq='D').values)
samplingpoints = list(df.SamplingPoint.unique())

new_idx = []
for sp in samplingpoints:
    for d in dates:
        new_idx.append((sp, np.datetime64(d)))

df.set_index(keys=['SamplingPoint', 'DatetimeBegin'], inplace=True)
df.sort_index(inplace=True)
df = df.reindex(new_idx)
#print(df.loc['SPO-BETR223_00001_100','2013-01-29'])  # should contain NaN for the columns

df['AirPollutionLevel'] = df.groupby(level=0).AirPollutionLevel.bfill().fillna(0)
#print(df.loc['SPO-BETR223_00001_100','2013-01-29'])  # NaN are replaced by values of 2013-01-30
print('{} missing values'.format(df.isnull().sum().sum()))

Handling multiple time series

Alright, now we have a data set that is cleaned and does not contain any missing values. One aspect that makes this data set particular is that we have data for multiple samplingpoints. So we have multiple time series.

One way to deal with that is to create dummy variables for the samplingpoints and use all records to train the model. Another way is to build a separate model per samplingpoint.

In this notebook, we will do the latter. We will, however, limit the notebook to do that for only one samplingpoint. But the same logic can be applied to every samplingpoint.

df = df.loc['SPO-BETR223_00001_100',:]

Split train, test and validation set

We split off a test set in order to evaluate the performance of the model. The test set will not be used during the training phase.

train set: data until July 2014
validation set: 6 months between July 2014 and January 2015
test set: data of 2015

train = df.query('DatetimeBegin < "2014-07-01"')
valid = df.query('DatetimeBegin >= "2014-07-01" and DatetimeBegin < "2015-01-01"')
test = df.query('DatetimeBegin >= "2015-01-01"')

Scaling

# Save column names and indices to use when storing as csv
cols = train.columns
train_idx = train.index
valid_idx = valid.index
test_idx = test.index

# normalize the dataset
scaler = MinMaxScaler(feature_range=(0, 1))
train = scaler.fit_transform(train)
valid = scaler.transform(valid)
test = scaler.transform(test)

Save the processed datasets

That way we don’t need to redo the preprocessing every time we rerun the notebook.

train = pd.DataFrame(train, columns=cols, index=train_idx)
valid = pd.DataFrame(valid, columns=cols, index=valid_idx)
test = pd.DataFrame(test, columns=cols, index=test_idx)

train.to_csv('../data/processed/train.csv')
valid.to_csv('../data/processed/valid.csv')
test.to_csv('../data/processed/test.csv')

Modeling

First, we read in the processed data sets. Secondly, we create a function to plot the training and validation loss for the different models we will build.

train = pd.read_csv('../data/processed/train.csv', header=0, index_col=0).values.astype('float32')
valid = pd.read_csv('../data/processed/valid.csv', header=0, index_col=0).values.astype('float32')
test = pd.read_csv('../data/processed/test.csv', header=0, index_col=0).values.astype('float32')

def plot_loss(history, title):
    plt.figure(figsize=(10,6))
    plt.plot(history.history['loss'], label='Train')
    plt.plot(history.history['val_loss'], label='Validation')
    plt.title(title)
    plt.xlabel('Nb Epochs')
    plt.ylabel('Loss')
    plt.legend()
    plt.show()

    val_loss = history.history['val_loss']
    min_idx = np.argmin(val_loss)
    min_val_loss = val_loss[min_idx]
    print('Minimum validation loss of {} reached at epoch {}'.format(min_val_loss, min_idx))

Prepare data with the TimeseriesGenerator

The TimeseriesGenerator of Keras helps us building the data in the correct format for modeling.

length: number of timesteps in the generated sequence. Here we want to look back an arbitrary number of _nlag timesteps. In reality, n_lag could depend on how the predictions will be used. Suppose the Belgian government can take some actions to reduce the SO2 pollution around a samplingpoint (for instance prohibit the entrance of diesel cars in a city for a certain period of time). And suppose the government needs 14 days before the corrective actions can go in effect. Then it would make sense to set n_lag to 14.
sampling_rate: number of timesteps between successive timesteps in the generated sequence. We want to keep all timesteps, so we leave this to the default of 1.
stride: this parameter influences how much the generated sequences will overlap. As we do not have much data, we leave it to the default of 1. This means that two sequences generated after one another overlap with all timesteps except one.
batch_size: number of generated sequences in each batch

n_lag = 14

train_data_gen = TimeseriesGenerator(train, train, length=n_lag, sampling_rate=1, stride=1, batch_size = 5)
valid_data_gen = TimeseriesGenerator(train, train, length=n_lag, sampling_rate=1, stride=1, batch_size = 1)
test_data_gen = TimeseriesGenerator(test, test, length=n_lag, sampling_rate=1, stride=1, batch_size = 1)

Recurrent Neural Networks

Traditional neural networks have no memory. Consequently, they do not take into account previous input when processing the current input. In sequential data sets, like time series, the information of previous time steps is typically relevant for predicting something in the current step. So a state about the previous time steps needs to be maintained.

In our case, the air pollution at time t might be influenced by air pollution in previous timesteps. So we need to take that into account. Recurrent Neural Networks or RNNs have an internal loop by which they maintain a state of previous timesteps. This state is then used for the prediction in the current timestep. The state is reset when a new sequence is being processed.

For an illustrated guide about RNNs, you should definitely read the article by Michael Nguyen.

For our case, we use a SimpleRNN of the Keras package. We also specify an EarlyStopping callback to stop the training when there were 10 epochs without any improvement in the validation loss. The ModelCheckpoint allows us to save the weights of the best model. The model architecture still needs to be saved separately.

simple_rnn = Sequential()
simple_rnn.add(SimpleRNN(4, input_shape=(n_lag, 1)))
simple_rnn.add(Dense(1))
simple_rnn.compile(loss='mae', optimizer=RMSprop())

checkpointer = ModelCheckpoint(filepath='../model/simple_rnn_weights.hdf5'
                               , verbose=0
                               , save_best_only=True)
earlystopper = EarlyStopping(monitor='val_loss'
                             , patience=10
                             , verbose=0)
with open("../model/simple_rnn.json", "w") as m:
    m.write(simple_rnn.to_json())

simple_rnn_history = simple_rnn.fit_generator(train_data_gen
                                              , epochs=100
                                              , validation_data=valid_data_gen
                                              , verbose=0
                                              , callbacks=[checkpointer, earlystopper])
plot_loss(simple_rnn_history, 'SimpleRNN - Train & Validation Loss')

Training and validation loss for a SimpleRNN

Long Short Term Memory Networks

An RNN has a short memory. It has difficulty remembering information from many timesteps ago. This occurs when the sequences are very long.

In fact, it is due to the vanishing gradient problem. The gradients are values that update the weights of a neural network. When you have many timesteps in your RNN the gradient for the first layers becomes very tiny. As a result, the update of the weights of the first layers is negligible. This means that the RNN is not capable of learning what was in the early layers.

So we need a way to carry the information of the first layers to later layers. LSTMs are better suited to take into account long-term dependencies. Michael Nguyen wrote an excellent article with a visual description of LSTMs.

Simple LSTM model

simple_lstm = Sequential()
simple_lstm.add(LSTM(4, input_shape=(n_lag, 1)))
simple_lstm.add(Dense(1))
simple_lstm.compile(loss='mae', optimizer=RMSprop())

checkpointer = ModelCheckpoint(filepath='../model/simple_lstm_weights.hdf5'
                               , verbose=0
                               , save_best_only=True)
earlystopper = EarlyStopping(monitor='val_loss'
                             , patience=10
                             , verbose=0)
with open("../model/simple_lstm.json", "w") as m:
    m.write(simple_lstm.to_json())

simple_lstm_history = simple_lstm.fit_generator(train_data_gen
                                                , epochs=100
                                                , validation_data=valid_data_gen
                                                , verbose=0
                                                , callbacks=[checkpointer, earlystopper])
plot_loss(simple_lstm_history, 'Simple LSTM - Train & Validation Loss')

Training and validation loss for a simple LSTM

Stacked LSTM model

In this model, we will be stacking multiple LSTM layers. That way the model will learn other abstractions of input data over time. In other words, representing the input data at different time scales.

To do that in Keras, we need to specify the parameter return_sequences in the LSTM layer preceding another LSTM layer.

stacked_lstm = Sequential()
stacked_lstm.add(LSTM(16, input_shape=(n_lag, 1), return_sequences=True))
stacked_lstm.add(LSTM(8, return_sequences=True))
stacked_lstm.add(LSTM(4))
stacked_lstm.add(Dense(1))
stacked_lstm.compile(loss='mae', optimizer=RMSprop())

checkpointer = ModelCheckpoint(filepath='../model/stacked_lstm_weights.hdf5'
                               , verbose=0
                               , save_best_only=True)
earlystopper = EarlyStopping(monitor='val_loss'
                             , patience=10
                             , verbose=0)
with open("../model/stacked_lstm.json", "w") as m:
    m.write(stacked_lstm.to_json())

stacked_lstm_history = stacked_lstm.fit_generator(train_data_gen
                                                  , epochs=100
                                                  , validation_data=valid_data_gen
                                                  , verbose=0
                                                  , callbacks=[checkpointer, earlystopper])
plot_loss(stacked_lstm_history, 'Stacked LSTM - Train & Validation Loss')

Training and validation loss for a stacked LSTM

Evaluating Performance

Based on the minimum validation losses the SimpleRNN seems to outperform the LSTM models, although the metrics are close to each other.

With the evaluate_generator method, we can evaluate the models on the test data (generator). This will give us the loss on the test data. We will first load the model architecture from the JSON files and the best model’s weights.

def eval_best_model(model):
    # Load model architecture from JSON
    model_architecture = open('../model/'+model+'.json', 'r')
    best_model = model_from_json(model_architecture.read())
    model_architecture.close()
    # Load best model's weights
    best_model.load_weights('../model/'+model+'_weights.hdf5')
    # Compile the best model
    best_model.compile(loss='mae', optimizer=RMSprop())
    # Evaluate on test data
    perf_best_model = best_model.evaluate_generator(test_data_gen)
    print('Loss on test data for {} : {}'.format(model, perf_best_model))

eval_best_model('simple_rnn')
eval_best_model('simple_lstm')
eval_best_model('stacked_lstm')

Loss on test data for simple_rnn: 0.01638169982337905
Loss on test data for simple_lstm: 0.015934137431135205
Loss on test data for stacked_lstm: 0.015420083056716116

Conclusion

In this story, we used a Recurrent Neural Network and two different architectures for an LSTM. The best performance comes from the stacked LSTM consisting of a few hidden layers.

There are definitely a number of things worth investigating further that could improve the model's performance.

Use the hourly data (another CSV file available on the EEA website) and try out other sampling strategies than the daily data.
Use the data about the other pollutants as features to predict SO2 pollution. Perhaps other pollutants are correlated to the SO2 pollution.
Construct other features based on the date. A nice write-up can be found in the PDF of one of the winners of the Power Laws Forecasting competition of Driven Data.

I’ve learned a lot about recurrent neural networks by doing this project. I hope you’ve enjoyed it. Feel free to leave any comments!