Precipitation Nowcasting se strojovým učením

V současné době meteorologové odhadují, že 90 % předpovědí počasí je správných během 5 dnů. Provedené předpovědi jsou obvykle založeny na dvou samostatných metodách:

Fyzikální přístupy: Tyto přístupy využívají modely, které integrují měřitelné veličiny, jako je tlak, pohyb mraků, podmínky na obloze... Takové modely jsou dobré při předpovídání počasí na nadcházející dny nebo týdny.
Fyzikální přístupy (založené na datech): Tyto přístupy využívají historická data k vytváření modelů, které mohou předpovídat. Takové modely vykazují dobré výsledky v předpovídání počasí až na 6 hodin nebo to, co je známé jako weather nowcasting.

V tomto článku se budeme zabývat druhou kategorií přístupů. Budeme diskutovat o různých formátech dat o počasí, o tom, jak lze využít strojové učení k provádění předpovědi počasí a jak mohou být vyvinuté architektury přínosné pro řešení podobných problémů.

Data pro předpověď klimatu

Protože přístupy bez fyziky používají historická data, začněme tím, že se podíváme na dostupná data.

Budeme používat dva hlavní zdroje dat:

Obrazová data: Tato data mají formu radarových nebo satelitních snímků určité zeměpisné oblasti. Slouží k předpovědi srážek, pohybu větru nebo vlhkosti.

Tabulková data : Tato data mají formu záznamů měřitelných veličin jako je teplota, vlhkost nebo rychlost větru.

Zatímco oba zdroje dat jsou důležité pro vytváření výkonných modelů, zaměříme se z důvodu jednoduchosti na první (obrázková data shromážděná z radarů nebo satelitů). Nejčastěji používané modely s obrazovými daty jsou Convolutional Neural Networks (CNN ).

Po této práci použijeme architekturu U-Net k vytvoření vlastního modelu nowcastingu.

Architektura

Vycházet ze stávající architektury je užitečné z mnoha důvodů, mezi které patří:

Architektury slouží jako vodítko pro vytváření nových modelů. Lidé, kteří vytvářejí nové architektury, používají metodu pokusů a omylů, aby dosáhli konečných výsledků. Opětovným použitím jejich konečných výsledků můžeme ušetřit spoustu času.

Předtrénované modely jsou obvykle k dispozici pro okamžité použití. Když výzkumníci publikují své nové architektury, obvykle zveřejňují také natrénované parametry, takže uživatelé nebudou muset procházet potížemi se školením/optimalizací od nuly. To je užitečné zejména pro velmi velké modely chtivé zdrojů.

Mezi příklady slavných architektur vidění patří:

LeNet (60k parametrů)

AlexNet (60m parametrů)

VGG-16 (138 m parametrů)

U-Net

U-Net je architektura založená na plně konvoluční síti, což znamená, že nemá žádné plně propojené vrstvy. Poprvé byl představen pro úlohu segmentace lékařského obrazu. Tyto výsledky inspirovaly výzkumníky k tomu, aby jej rozšířili na další úkoly v počítačovém vidění.

V roce 2019 použil Google architekturu založenou na U-Net k vytvoření modelu předpovědi srážek.

Název „U-Net“ pochází z tvaru „U“ jeho architektury.

Rozlišujeme tři hlavní složky:

Contracting / Encoder: Posloupnost vrstev konvoluce/sdružování, které komprimují vstupní obraz do reprezentace menší velikosti.
Most/úzké hrdlo: Spodní část „U“ spojuje kodér s dekodérem. Je tvořena řadou konvolučních operací.
Dekontrakce / dekodér : Posloupnost vzestupných a konvolučních vrstev, které „dekomprimují“ výstup úzkého hrdla.

Architektura U-Net zatím vypadá jako automatický kodér. Rozdíl však spočívá v předávání informací mezi kodérem a dekodérem. Tyto informace jsou předávány zřetězením výsledků konvolucí z kodéru s výsledky vzestupné konvoluce dekodéru. Tato úprava zvyšuje rozlišení a umožňuje modelu vydávat prostorově přesnější výstup (umístění pixelů ve výstupu bude přesnější, což řeší problém, který se vyskytl u modelů, které neobsahovaly zřetězení). Zřetězení se provádí symetricky.

Vytvoření modelu U-Net pro srážkový Nowcasting

V této části vytvoříme model U-Net pro předpovídání počasí ze satelitních snímků. Použijeme předtrénovaný model nazvaný Rain-Net.

Níže uvedený kód je k dispozici v tomto colab.

Nejprve nainstalujeme wradlib, knihovnu s otevřeným zdrojovým kódem pro zpracování dat meteorologických radarů

!pip install wradlib
import wradlib as wrl

Poté napíšeme dvě pomocné funkce pro stažení satelitních dat z otevřeného datového serveru 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

Tyto nástroje nám umožňují stáhnout nejnovější 4 satelitní snímky, což je počet satelitních snímků, které potřebujeme k předpovědi pomocí našeho předem připraveného modelu.

Vytvořené funkce pak můžeme použít k získání nejnovějších 4 obrázků


RY_latest, RY_latest_timestep = download_data()

Po získání obrázků použijeme metodu wradlib vis.plot_ppi k vykreslení dat

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

Nyní jsme načetli naše data. Dále načteme model.

Začneme importem příslušných tříd. V tomto článku budeme používat TensorFlow.

z tensorflow.keras.layers import Input, Conv2D, Activation, Concatenate, Conv2DTransose, MaxPool2D

z tensorflow.keras.models import Model

Postavme 3 primitivní stavební bloky. Tyto „stavební bloky“ budou použity k vytvoření celé architektury podle této implementace.

První blok odpovídá posloupnosti konvolučních vrstev, nazýváme ho "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

Druhý blok se používá k sestavení části kodéru (konvoluční blok + maximální sdružování). Říkáme tomu „kódovací_blok“

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

Třetí a poslední blok je „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

Kombinujeme tyto stavební bloky a vytváříme model 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

Klidně si s implementací pohrajte, vylepšete ji nebo si ji přizpůsobte svým potřebám, pokud ji chcete použít na něco jiného.

Trénink tohoto modelu může zabrat hodně času. Naštěstí existuje model nazvaný Rain-Net, který byl vytvořen na základě architektury U-Net a specializuje se na srážky nowcasting.

Pojďme naklonovat jeho úložiště GitHub

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

Následně si stáhneme předtrénované závaží pro tento model

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

Dalším krokem je vytvoření modelu založeného na architektuře nalezené v úložišti a poté načtení vah stažených do tohoto modelu

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

Obrázky, které jsme stáhli, mají velikost 900*900 pixelů. Chystáme se přetvořit tyto obrázky tak, aby odpovídaly očekávanému vstupu 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

Poté vytvoříme funkci, která předpovídá.

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

Tuto funkci pak voláme na datech, která jsme stáhli dříve

Y_pred = prediction(model, RY_latest)

Můžeme vykreslit předpovědi a uložit je, abychom z uložených výsledků vytvořili obrázek gif, který nám umožní vizualizovat předpovědi

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

Gratulujeme, že jste to dotáhli tak daleko! Nyní můžete používat Rain-Net k vytváření předpovědí a jejich vizualizaci.

Závěr

V tomto článku jsme použili model strojového učení (Rain-Net) k provádění srážek nowcasting. Použili jsme architekturu U-Net, kterou jsme postavili pomocí TensorFlow. Nahráli jsme předtrénovaný model pro předpovídání satelitních snímků.

Tuto implementaci lze zlepšit mnoha způsoby. Například:

Dolaďte model na datové sadě
V architektuře použijte modul založený na pozornosti, jako je CBAM (Convolutional Block Attention Module).

Reference

– Colab

Původní papír U-Net

– Blog Googlu o jejich modelu předpovědi počasí založeném na U-Net

Implementace Rain-Net

– Implementace U-Net

Přijďte na jeden z našich bezplatných workshopů!

Začněte svou kariéru datového vědce pomocí našich bezplatných workshopů, které jsou založeny na přizpůsobitelném učebním plánu a vedené odborníky z oboru.