Precipitații Nowcasting cu Machine Learning

În zilele noastre, meteorologii estimează că 90% din prognozele meteo sunt corecte într-un interval de 5 zile. Predicțiile care sunt făcute se bazează de obicei pe două metode separate:

Abordări bazate pe fizică: Aceste abordări folosesc modele care integrează cantități măsurabile, cum ar fi presiunea, mișcarea norilor, condițiile cerului... Astfel de modele sunt bune la prezicerea vremii pentru zilele sau săptămânile următoare.
Abordări fără fizică (bazate pe date): Aceste abordări folosesc date istorice pentru a crea modele care pot face predicții. Astfel de modele arată rezultate bune în prezicerea vremii pentru până la 6 ore sau ceea ce este cunoscut sub numele de prognoza meteo.

În acest articol, vom discuta despre a doua categorie de abordări. Vom discuta despre diferitele formate de date meteorologice, despre modul în care învățarea automată poate fi valorificată pentru a face prognoza meteo și despre modul în care arhitecturile dezvoltate pot fi benefice pentru rezolvarea unor probleme similare.

Date pentru predicția climei

Deoarece abordările fără fizică folosesc date istorice, să începem prin a analiza datele disponibile.

Vom folosi două surse principale de date:

Date de imagine: Aceste date iau forma unor imagini radar sau satelit ale unei zone geografice specificate. Este folosit pentru a prezice precipitații, mișcarea vântului sau umiditatea.

Date tabelare: Aceste date iau forma unor înregistrări ale cantităților măsurabile, cum ar fi temperatura, umiditatea sau viteza vântului.

Deși ambele surse de date sunt importante pentru a construi modele puternice, ne vom concentra pe primele (date de imagine colectate de la radare sau sateliți) din motive de simplitate. Cele mai frecvent utilizate modele cu date de imagine sunt Convolutional Neural Networks (CNNs ).

În urma acestei lucrări, vom folosi o arhitectură U-Net pentru a ne construi propriul model nowcasting.

Arhitectură

Pornirea de la o arhitectură existentă este utilă din mai multe motive, printre care:

Arhitecturile servesc drept linii directoare pentru crearea de noi modele. Oamenii care creează noi arhitecturi adoptă o abordare de încercare și eroare pentru a ajunge la rezultatele finale. Prin reutilizarea rezultatelor lor finale, putem economisi mult timp.

Modelele pregătite sunt de obicei disponibile pentru utilizare imediată. Când cercetătorii își publică noile arhitecturi, de obicei publică și parametrii antrenați, astfel încât utilizatorii să nu fie nevoiți să treacă prin necazul antrenării/optimizării de la zero. Acest lucru este util în special pentru modele foarte mari, însetate de resurse.

Exemple de arhitecturi de viziune celebre includ:

LeNet (60k parametri)

AlexNet (parametri 60m)

VGG-16 (parametri 138m)

U-Net

U-Net este o arhitectură bazată pe o rețea complet convoluțională, ceea ce înseamnă că nu are straturi complet conectate. A fost introdus pentru prima dată pentru o sarcină de segmentare a imaginii medicale. Aceste rezultate i-au inspirat pe cercetători să-l extindă la alte sarcini din viziunea computerizată.

În 2019, Google a folosit o arhitectură bazată pe U-Net pentru a crea un model de prognoză a precipitațiilor.

Numele „U-Net” provine de la forma „U” a arhitecturii sale.

Identificăm trei componente principale:

Contractare / Encoder: Succesiune de straturi de convoluție/pooling care comprimă imaginea de intrare într-o reprezentare de dimensiune mai mică.
Pod/Gât: Partea de jos a „U” conectează codificatorul la decodor. Este format dintr-o serie de operații de convoluție.
Decontractare / Decodor: Succesiunea de convoluții în sus și straturi de convoluție, care „decomprimă” ieșirea blocajului.

Arhitectura unui U-Net arată ca un auto-encoder până acum. Totuși, diferența constă în informațiile care trec între codificator și decodor. Această informație este transmisă prin concatenarea rezultatelor convoluțiilor de la codificator cu cele ale convoluției în sus a decodorului. Această modificare îmbunătățește rezoluția și permite modelului să scoată o ieșire mai precisă din punct de vedere spațial (locația pixelilor din ieșire va fi mai precisă, rezolvând o problemă care a apărut în modelele care nu includeau concatenare). Concatenarea se face simetric.

Construirea unui model U-Net pentru Precipitații Nowcasting

În această secțiune, vom construi un model U-Net pentru a prezice vremea din imagini din satelit. Vom folosi un model pre-antrenat numit Rain-Net.

Codul de mai jos este disponibil în acest colab.

Mai întâi instalăm wradlib, o bibliotecă Open Source pentru procesarea datelor radar meteorologic

!pip install wradlib
import wradlib as wrl

Apoi scriem două funcții utilitare pentru a descărca datele satelitului de pe serverul de date deschise al 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

Aceste utilități ne permit să descarcăm cele mai recente 4 imagini din satelit, care este numărul de imagini din satelit de care avem nevoie pentru a face predicții cu modelul nostru pregătit.

Putem folosi apoi funcțiile create pentru a obține cele mai recente 4 imagini


RY_latest, RY_latest_timestep = download_data()

După obținerea imaginilor, folosim metoda wradlib vis.plot_ppi pentru a reprezenta datele

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

Acum ne-am încărcat datele. Să încărcăm modelul în continuare.

Începem prin a importa clasele relevante. Vom folosi TensorFlow în acest articol.

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

din tensorflow.keras.models import Model

Să construim 3 blocuri de construcție primitive. Aceste „blocuri” vor fi folosite pentru a crea întreaga arhitectură, conform acestei implementare.

Primul bloc corespunde succesiunii straturilor convoluționale, îl numim „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

Al doilea bloc este folosit pentru a construi partea din codificator (bloc convoluțional + pooling maxim). Îl numim „coder_block”

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

Al treilea și ultimul bloc este un „decoder_block” (convoluție în sus + concatenare + convoluție).

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

Combinăm aceste blocuri pentru a construi modelul 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

Simțiți-vă liber să vă jucați cu implementarea, să o îmbunătățiți sau să o adaptați nevoilor dvs. dacă doriți să o utilizați pentru altceva.

Antrenarea acestui model poate dura mult timp. Din fericire, există un model numit Rain-Net, care a fost creat pe baza unei arhitecturi U-Net și este specializat în prognozarea precipitațiilor.

Să clonăm depozitul său GitHub

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

Descărcăm apoi greutățile preantrenate pentru acest model

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

Următorul pas este să creați un model bazat pe arhitectura găsită în depozit, apoi să încărcați greutățile descărcate în acest model.

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

Imaginile pe care le-am descărcat au o dimensiune de 900*900 pixeli. Vom remodela aceste imagini pentru a se potrivi cu intrarea așteptată a 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

Apoi creăm o funcție care face predicțiile.

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

Apoi apelăm această funcție pe datele pe care le-am descărcat mai devreme

Y_pred = prediction(model, RY_latest)

Putem trasa predicțiile și le salva pentru a folosi rezultatele salvate pentru a crea o imagine gif care ne permite să vizualizăm predicțiile

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

Felicitări pentru că ai ajuns până aici! Acum puteți folosi un Rain-Net pentru a face predicții și a le vizualiza.

Concluzie

În acest articol, am folosit un model de învățare automată (Rain-Net) pentru a face previzionarea precipitațiilor. Am folosit o arhitectură U-Net, pe care am construit-o folosind TensorFlow. Am încărcat un model preantrenat pentru a prezice imagini din satelit.

Această implementare poate fi îmbunătățită în multe feluri. De exemplu:

Ajustați modelul pe setul dvs. de date
Utilizați un modul bazat pe atenție, cum ar fi CBAM (Convolutional Block Attention Module) în arhitectură.

Referințe

Vino la unul dintre atelierele noastre gratuite!

Începeți-vă cariera ca cercetător de date cu atelierele noastre gratuite, care se bazează pe un curriculum adaptabil și sunt ghidate de experți din industrie.