Nowcasting sulle precipitazioni con il machine learning

Al giorno d'oggi, i meteorologi stimano che il 90% delle previsioni meteorologiche siano corrette in un arco di 5 giorni. Le previsioni che vengono fatte si basano solitamente su due metodi distinti:

Approcci basati sulla fisica: questi approcci utilizzano modelli che integrano quantità misurabili come pressione, movimento delle nuvole, condizioni del cielo... Tali modelli sono efficaci nel prevedere il tempo per i prossimi giorni o settimane.
Approcci privi di fisica (basati sui dati): questi approcci utilizzano dati storici per creare modelli in grado di fare previsioni. Tali modelli mostrano buoni risultati nella previsione del tempo fino a 6 ore o nel cosiddetto nowcasting meteorologico.

In questo articolo discuteremo della seconda categoria di approcci. Discuteremo dei diversi formati dei dati meteorologici, di come sfruttare l'apprendimento automatico per effettuare previsioni meteorologiche e di come le architetture sviluppate possano essere utili per risolvere problemi simili.

Dati per la previsione climatica

Poiché gli approcci privi di fisica utilizzano dati storici, iniziamo esaminando i dati disponibili.

Utilizzeremo due principali fonti di dati:

Dati immagine: questi dati assumono la forma di immagini radar o satellitari di un'area geografica specifica. Viene utilizzato per prevedere le precipitazioni, il movimento del vento o l'umidità.

Dati tabulari: questi dati assumono la forma di registrazioni di quantità misurabili come temperatura, umidità o velocità del vento.

Sebbene entrambe le fonti di dati siano importanti per costruire modelli potenti, ci concentreremo sulla prima (dati di immagini raccolti da radar o satelliti) per ragioni di semplicità. I modelli più comunemente utilizzati con dati di immagine sono Reti neurali convoluzionali (CNN ).

A seguito di questo lavoro, utilizzeremo un'architettura U-Net per costruire il nostro modello di nowcasting.

Architettura

Partire da un’architettura esistente è utile per molteplici ragioni, tra cui:

Le architetture servono come linee guida per la creazione di nuovi modelli. Le persone che creano nuove architetture adottano un approccio per tentativi ed errori per arrivare ai risultati finali. Riutilizzando i risultati finali, possiamo risparmiare molto tempo.

I modelli preaddestrati sono generalmente disponibili per l'uso immediato. Quando i ricercatori pubblicano le loro nuove architetture, di solito pubblicano anche i parametri addestrati, in modo che gli utenti non debbano affrontare il problema dell'addestramento/ottimizzazione da zero. Ciò è particolarmente utile per modelli molto grandi e assetati di risorse.

Esempi di architetture di visione famose includono:

LeNet (parametri 60k)

AlexNet (parametri 60m)

VGG-16 (parametri 138m)

##U-Net

U-Net è un'architettura basata su una rete completamente convoluzionale, il che significa che non ha livelli completamente connessi. È stato introdotto per la prima volta per un'attività di segmentazione di immagini mediche. Questi risultati hanno ispirato i ricercatori ad estenderlo ad altri compiti della visione artificiale.

Nel 2019, Google ha utilizzato un'architettura basata su U-Net per creare un modello di previsione delle precipitazioni.

Il nome "U-Net" deriva dalla forma a "U" della sua architettura.

Identifichiamo tre componenti principali:

Contrazione/Codificatore: successione di strati di convoluzione/pooling che comprimono l'immagine di input in una rappresentazione di dimensioni più piccole.
Ponte/Collo di bottiglia: la parte inferiore della "U" collega l'encoder al decoder. È formato da una serie di operazioni di convoluzione.
Decontratturante/Decodificatore: successione di upconvoluzioni e strati di convoluzione, che “decomprimono” l'output del collo di bottiglia.

L'architettura di una U-Net finora sembra un codificatore automatico. Tuttavia, la differenza risiede nelle informazioni che passano tra il codificatore e il decodificatore. Queste informazioni vengono passate concatenando i risultati delle convoluzioni del codificatore con quelli della convoluzione superiore del decodificatore. Questa modifica migliora la risoluzione e consente al modello di produrre un output più preciso dal punto di vista spaziale (la posizione dei pixel nell'output sarà più precisa, risolvendo un problema che si verificava nei modelli che non includevano la concatenazione). La concatenazione avviene simmetricamente.

Costruire un modello U-Net per il nowcasting delle precipitazioni

In questa sezione costruiremo un modello U-Net per prevedere il tempo dalle immagini satellitari. Utilizzeremo un modello pre-addestrato chiamato Rain-Net.

Il codice seguente è disponibile in questo colab.

Per prima cosa installiamo wradlib, una libreria open source per l'elaborazione dei dati radar meteorologici

!pip install wradlib
import wradlib as wrl

Scriviamo quindi due funzioni di utilità per scaricare i dati satellitari dal server dati aperto di DWD


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

Queste utilità ci consentono di scaricare le ultime 4 immagini satellitari, ovvero il numero di immagini satellitari di cui abbiamo bisogno per fare previsioni con il nostro modello pre-addestrato.

Possiamo quindi utilizzare le funzioni create per ottenere le ultime 4 immagini


RY_latest, RY_latest_timestep = download_data()

Dopo aver ottenuto le immagini, utilizziamo il metodo vis.plot_ppi di wradlib per tracciare i dati

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

Ora abbiamo caricato i nostri dati. Carichiamo il modello successivo.

Iniziamo importando le classi rilevanti. Utilizzeremo TensorFlow in questo articolo.

da tensorflow.keras.layers importa Input, Conv2D, Attivazione, Concatena, Conv2DTranspose, MaxPool2D

da tensorflow.keras.models modello di importazione

Costruiamo 3 elementi costitutivi primitivi. Questi "mattoni" verranno utilizzati per creare l'intera architettura, secondo questa implementazione.

Il primo blocco corrisponde alla successione degli strati convoluzionali, lo chiamiamo "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

Il secondo blocco viene utilizzato per costruire dalla parte codificatore (blocco convoluzionale + pooling massimo). Lo chiamiamo “encoder_block”

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

Il terzo e ultimo blocco è un “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

Combiniamo questi elementi costitutivi per costruire il modello U-Net

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

Sentiti libero di giocare con l'implementazione, migliorarla o adattarla alle tue esigenze se vuoi usarla per qualcos'altro.

L'addestramento di questo modello può richiedere molto tempo. Fortunatamente, esiste un modello chiamato Rain-Net, che è stato creato sulla base di un’architettura U-Net ed è specializzato nel nowcasting delle precipitazioni.

Cloniamo il suo repository GitHub

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

Quindi scarichiamo i pesi preallenati per questo modello

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

Il passo successivo è creare un modello basato sull'architettura trovata nel repository, quindi caricare i pesi scaricati su questo modello

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

Le immagini che abbiamo scaricato hanno una dimensione di 900*900 pixel. Rimodelleremo queste immagini per adattarle all'input previsto di 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

Creiamo quindi una funzione che effettua le previsioni.

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

Chiamiamo quindi questa funzione sui dati che abbiamo scaricato in precedenza

Y_pred = prediction(model, RY_latest)

Possiamo tracciare i pronostici e salvarli per utilizzare i risultati salvati per creare un'immagine gif che ci permetta di visualizzare i pronostici

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')

Congratulazioni per essere arrivato fin qui! Ora puoi utilizzare una Rain-Net per fare previsioni e visualizzarle.

Conclusione

In questo articolo, abbiamo utilizzato un modello di apprendimento automatico (Rain-Net) per effettuare il nowcasting delle precipitazioni. Abbiamo utilizzato un'architettura U-Net, che abbiamo creato utilizzando TensorFlow. Abbiamo caricato un modello preaddestrato per prevedere le immagini satellitari.

Questa implementazione può essere migliorata in molti modi. Per esempio:

Ottimizza il modello sul tuo set di dati
Utilizzare un modulo basato sull'attenzione, come CBAM (Modulo di attenzione del blocco convoluzionale) nell'architettura.

Riferimenti

Vieni a uno dei nostri workshop gratuiti!

Inizia la tua carriera come data scientist con i nostri workshop gratuiti, basati su un curriculum adattabile e guidati da esperti del settore.