Precipitation Nowcasting with Machine Learning

Nowadays, meteorologists estimate that 90% of the weather predictions are correct in a 5 day span. The predictions that are made are usually based on two separate methods :

Physics-based approaches: These approaches use models that integrate measurable quantities such as pressure, clouds’ movement, sky conditions… Such models are good at predicting the weather for the upcoming days or weeks.
Physics-free (data-based) approaches: These approaches use historical data to create models that can make predictions. Such models show good results in predicting the weather for up to 6 hours or what is known as weather nowcasting.

In this article, we will discuss the second category of approaches. We will be discussing the different formats of weather data, how machine learning can be leveraged to do weather nowcasting, and how the architectures developed can be beneficial for solving similar problems.

Data for Climate Prediction

Since the physics-free approaches use historical data, let’s start by looking into the data available.

We will be using two main sources of data:

Image data: This data takes the form of radar or satellite images of a specified geographic area. It is used to predict precipitation, wind movement or humidity.

Tabular data : This data takes the form of records of measurable quantities such as temperature, humidity or wind speed.

While both data sources are important to build powerful models, we will be focusing on the former (image data collected from radars or satellites) for simplicity reasons. The most commonly used models with image data are Convolutional Neural Networks (CNNs).

Following this work, we are going to use a U-Net architecture to build our own nowcasting model.

Architecture

Starting from an existing architecture is helpful for many reasons, among which:

Architectures serve as guidelines for creating new models. People who create new architectures adopt a trial and error approach to get to their final results. By reusing their final results, we can save a lot of time.

Pretrained models are usually available for immediate use. When researchers publish their new architectures, they usually publish the trained parameters as well, so that users won’t have to go through the trouble of training/optimizing from scratch. This is especially useful for very large, resource-thirsty models.

Example famous vision architectures include:

LeNet (60k parameters)

AlexNet (60m parameters)

VGG-16 (138m parameters)

U-Net

U-Net is an architecture based on a fully convolutional network, meaning it doesn’t have any fully connected layers. It was first introduced for a medical image segmentation task. These results inspired researchers to extend it to other tasks in computer vision.

In 2019, Google used a U-Net-based architecture to create a precipitation forecasting model.

The name “U-Net” comes from the “U” shape of its architecture.

We identify three main components:

Contracting / Encoder: Succession of convolution/pooling layers which compress the input image into a smaller size representation.
Bridge / Bottleneck: The bottom part of the “U” connects the encoder to the decoder. It is formed by a series of convolution operations.
Decontracting / Decoder : Succession of upconvolutions and convolution layers, which “decompress” the output of the bottleneck.

The architecture of a U-Net looks like an auto-encoder so far. However, the difference resides in the information passing between the encoder and the decoder. This information is passed by concatenating the results of convolutions from the encoder with that of the upconvolution of the decoder. This modification enhances the resolution, and allows for the model to output a more spatially precise output (the location of the pixels in the output will be more precise, solving a problem that occurred in models that didn’t include concatenation). The concatenation is done symmetrically.

Building a U-Net Model for Precipitation Nowcasting

In this section, we will build a U-Net model to predict the weather from satellite images. We will use a pre-trained model called Rain-Net.

The code below is available in this colab.

We first install wradlib, an Open Source Library for Weather Radar Data Processing

!pip install wradlib
import wradlib as wrl

We then write two utility functions to download satellite data from the DWD’s open data server


import urllib.request
import io
import numpy as np
import datetime

def download_and_read_RY(RY_timestamp):
   url = f"https://opendata.dwd.de/weather/radar/radolan/ry/raa01-ry_10000-{RY_timestamp}-dwd---bin"
   data_binary = urllib.request.urlopen(url).read()
   data, attr = wrl.io.read_radolan_composite(io.BytesIO(data_binary), missing = 0)
   data = data.astype("float32")
   return data,attr

def download_data():
   latest_scan, latest_attr = download_and_read_RY("latest")
   latest_datetime = latest_attr["datetime"]
   list_for_downloading = [ts.strftime("%y%m%d%H%M") for ts in
       [latest_datetime - datetime.timedelta(minutes=t) for t in [15, 10, 5]]]
   previous_scans = np.array([download_and_read_RY(ts)[0] for ts in list_for_downloading])
   print(list_for_downloading)
   scans = np.concatenate([previous_scans, latest_scan[np.newaxis, ::, ::]], axis = 0)

   return scans, latest_datetime

These utilities allow us to download the latest 4 satellite images, which is the number of satellite images that we need to make predictions with our pretrained model.

We can then use the functions created to get the latest 4 images


RY_latest, RY_latest_timestep = download_data()

After getting the images, we use wradlib’s vis.plot_ppi method to plot the data

for i in range(RY_latest.shape[0]):
   wrl.vis.plot_ppi(RY_latest[i])

We now have loaded our data. Let’s load the model next.

We start by importing the relevant classes. We will be using TensorFlow in this article.

from tensorflow.keras.layers import Input, Conv2D, Activation, Concatenate, Conv2DTranspose, MaxPool2D

from tensorflow.keras.models import Model

Let’s construct 3 primitive building blocks. These "building blocks" will be used to create the whole architecture, per this implementation.

The first block corresponds to the succession of convolutional layers, we call it "conv_block"

def conv_block(input, num_filters):
   x = Conv2D(num_filters, 3, padding ="same", kernel_initializer="he_normal"(input))
   x = Activation("relu")(x)
   x = Conv2D(num_filters, 3, padding ="same", kernel_initializer="he_normal"(x))
   x = Activation("relu")(x)
   return x

The second block is used to build the from the encoder part (convolutional block + max pooling). We call it an “encoder_block”

def encoder_block(input, num_filters):
   x = conv_block(input, num_filters)
   p = MaxPool2D((2,2))(x)
   return x,p

The third and final block is a “decoder_block” (upconvolution + concatenation + convolution).

def decoder_block(input, skip_features, num_filters):
   x = Conv2DTranspose(num_filters, (2,2), strides=2, padding="same", kernel_initializer="he_normal")
   x = Concatenate()([x, skip_features])
   x = conv_block(x, num_filters)
   return x

We combine these building blocks to construct the U-Net model

def build_unet(input_shape = (928, 928, 4)):

   inputs = Input(input_shape)
   e1, p1 = encoder_bock(inputs, 64)
   e2, p2 = encoder_bock(p1, 128)
   e3, p3 = encoder_bock(p2, 256)
   e4, p4 = encoder_bock(p3, 512)

   b1 = conv_block(p4, 1024)

   d1 = decoder_block(b1, e4, 512)
   d2 = decoder_block(d1, e3, 256)
   d3 = decoder_block(d2, e2, 128)
   d4 = decoder_block(d3, e1, 64)

   outputs = Conv2D(1, 1, padding="same", activation="linear")(d4)

   unet = Model(inputs, outputs, name="U-Net")

   return unet

Feel free to play around with the implementation, improve it or adapt it to your needs if you want to use it for something else.

Training this model can take a lot of time. Luckily, there is a model called Rain-Net, that has been created based on a U-Net architecture, and is specialized in precipitation nowcasting.

Let’s clone its GitHub repository

! git clone https://github.com/hydrogo/rainnet.git

We then download the pretrained weights for this model

!wget -O /content/rainnet/rainnet_weights.h5 https://zenodo.org/record/3630429/files/rainnet_weights.h5

The next step is to create a model based on the architecture found on the repository then load the weights downloaded to this model

import sys
from rainnet import rainnet
model = rainnet.rainnet()
model.load_weights('/content/rainnet/rainnet_weights.h5')
model.summary()

The images we downloaded have a size of 900*900 pixels. We are going to reshape these images to match the expected input of Rain-Net

def Scaler(array):
   return np.log(array+0.01)


def invScaler(array):
   return np.exp(array) - 0.01


def pad_to_shape(array, from_shape=900, to_shape=928, how="mirror"):
   # calculate how much to pad in respect with native resolution
   padding = int( (to_shape - from_shape) / 2)
   # for input shape as (batch, W, H, channels)
   if how == "zero":
       array_padded = np.pad(array, ((0,0),(padding,padding),(padding,padding),(0,0)), mode="constant", constant_values=0)
   elif how == "mirror":
       array_padded = np.pad(array, ((0,0),(padding,padding),(padding,padding),(0,0)), mode="reflect")
   return array_padded


def pred_to_rad(pred, from_shape=928, to_shape=900):
   # pred shape 12,928,928
   padding = int( (from_shape - to_shape) / 2)
   return pred[::, padding:padding+to_shape, padding:padding+to_shape].copy()

'''
the spatial extent of input data has to be a multiple of 2n+1
where n is the number of max pooling layers
'''

def data_preprocessing(X):

   # 0. Right shape for batch
   X = np.moveaxis(X, 0, -1)
   X = X[np.newaxis, ::, ::, ::]
   # 1. To log scale
   '''
   We use a log scale to respond to skewness towards large values
   '''
   X = Scaler(X)
   # 2. from 900x900 to 928x928
   X = pad_to_shape(X)

   return X


def data_postprocessing(nwcst):

   # 0. Squeeze empty dimensions
   nwcst = np.squeeze(np.array(nwcst))
   # 1. Convert back to rainfall depth
   nwcst = invScaler(nwcst)
   # 2. Convert from 928x928 back to 900x900
   nwcst = pred_to_rad(nwcst)
   # 3. Return only positive values
   nwcst = np.where(nwcst>0, nwcst, 0)
   return nwcst

We then create a function that makes the predictions.

def prediction(model_instance, input_data, lead_time=24):

   input_data = data_preprocessing(input_data)
   nwcst = []

   for _ in range(lead_time):
       pred = model_instance.predict(input_data)
       nwcst.append(pred)
       input_data = np.concatenate([input_data[::, ::, ::, 1:], pred], axis = -1)

   nwcst = data_postprocessing(nwcst)

   return nwcst

We then call this function on the data we downloaded earlier

Y_pred = prediction(model, RY_latest)

We can plot the predictions, and save them to use the results saved to create a gif image that allows us to visualize the predictions

import matplotlib.pyplot as plt
names = []
for i in range(0, 19):
   im = wrl.vis.plot_ppi(Y_pred[i])
   name = '/content/image'+str(i)+'.png'
   names.append(name)
   plt.savefig(name)

import imageio
from google.colab import files

imgs = []

for name in names:
   imgs.append(imageio.imread(name))

imageio.mimsave('/content/pred.gif', imgs, fps=4)
files.download('/content/pred.gif')

Congratulations for making it this far! You are now able to use a Rain-Net to make predictions and visualize them.

Conclusion

In this article, we have used a machine learning model (Rain-Net) to do precipitation nowcasting. We used a U-Net architecture, which we built using TensorFlow. We loaded a pretrained model to predict satellite images.

This implementation can be improved in many ways. For example:

Finetune the model on your dataset
Use an attention-based module, such as CBAM (Convolutional Block Attention Module) in the architecture.

References

Come To One Of Our Free Workshops!

Start your career as a data scientist with our free workshops, which are based on an adaptable curriculum and guided by industry experts.