機械学習による降水量ナウキャスティング

現在、気象学者は、5 日間の天気予報の 90% が正しいと推定しています。通常、行われる予測は 2 つの別個の方法に基づいています。

物理ベースのアプローチ: これらのアプローチでは、圧力、雲の動き、空の状態などの測定可能な量を統合するモデルが使用されます。このようなモデルは、今後数日または数週間の天気を予測するのに適しています。
物理学を使用しない (データベースの) アプローチ: これらのアプローチは、履歴データを使用して、予測を行うことができるモデルを作成します。このようなモデルは、最大 6 時間の天気を予測したり、いわゆる天気ナウキャスティングで優れた結果を示します。

この記事では、2 番目のカテゴリのアプローチについて説明します。さまざまな形式の気象データ、機械学習を活用して天気ナウキャスティングを行う方法、開発されたアーキテクチャが同様の問題の解決にどのように役立つかについて説明します。

気候予測のためのデータ

物理学を使用しないアプローチでは過去のデータが使用されるため、利用可能なデータを調べることから始めましょう。

次の 2 つの主要なデータソースを使用します。

画像データ: このデータは、指定された地理的領域のレーダーまたは衛星画像の形式をとります。降水量、風の動き、湿度を予測するために使用されます。

表形式データ : このデータは、温度、湿度、風速などの測定可能な量の記録の形式をとります。

強力なモデルを構築するには両方のデータソースが重要ですが、簡単にするために前者 (レーダーまたは衛星から収集された画像データ) に焦点を当てます。画像データで最も一般的に使用されるモデルは、畳み込みニューラルネットワーク (CNN) です。）。

この作業に従って、U-Net アーキテクチャを使用して独自のナウキャスティングモデルを構築します。

＃＃建築

既存のアーキテクチャから始めることは、次のような多くの理由で役立ちます。

アーキテクチャは、新しいモデルを作成するためのガイドラインとして機能します。新しいアーキテクチャを作成する人は、最終結果を得るために試行錯誤のアプローチを採用します。最終結果を再利用することで、時間を大幅に節約できます。

事前トレーニングされたモデルは通常、すぐに使用できます。研究者が新しいアーキテクチャを公開するときは、通常、トレーニング済みのパラメータも公開するため、ユーザーは最初からトレーニング/最適化を行う手間を省くことができます。これは、非常に大規模でリソースを大量に消費するモデルに特に役立ちます。

有名なビジョンアーキテクチャの例には次のものがあります。

LeNet (60k パラメータ)

AlexNet (60mパラメータ)

VGG-16 (138mパラメータ)

ユーネット

U-Net は完全畳み込みネットワークに基づくアーキテクチャです。つまり、完全に接続された層がありません。これは、医療画像のセグメンテーションタスクのために初めて導入されました。これらの結果に触発されて、研究者はそれをコンピュータービジョンの他のタスクにも拡張するようになりました。

2019 年、Google は U-Net ベースアーキテクチャを使用して降水量予測モデルを作成しました。

「U-Net」という名前は、そのアーキテクチャの「U」字型に由来しています。

3 つの主要なコンポーネントを特定します。

縮小/エンコーダー: 入力画像をより小さいサイズの表現に圧縮する畳み込み/プーリング層の連続。
ブリッジ/ボトルネック: 「U」の下部はエンコーダーとデコーダーを接続します。一連の畳み込み演算によって形成されます。
デコントラクト / デコーダー : ボトルネックの出力を「解凍」する、アップコンボリューションとコンボリューション層の連続。

U-Net のアーキテクチャは、これまでのところ自動エンコーダのように見えます。ただし、違いはエンコーダとデコーダの間で受け渡される情報にあります。この情報は、エンコーダからの畳み込みの結果とデコーダのアップコンボリューションの結果を連結することによって渡されます。この変更により解像度が向上し、モデルが空間的により正確な出力を出力できるようになります (出力内のピクセルの位置がより正確になり、連結を含まないモデルで発生した問題が解決されます)。連結は対称的に行われます。

降水量ナウキャスト用の U-Net モデルの構築

このセクションでは、衛星画像から天気を予測する U-Net モデルを構築します。 Rain-Net という事前トレーニング済みモデルを使用します。

以下のコードは、この colab で入手できます。

まず、気象レーダーデータ処理用のオープンソースライブラリである wradlib をインストールします。

!pip install wradlib
import wradlib as wrl

次に、DWD のオープンデータサーバーから衛星データをダウンロードする 2 つのユーティリティ関数を作成します。


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

これらのユーティリティを使用すると、最新の 4 つの衛星画像をダウンロードできます。これは、事前トレーニングされたモデルで予測を行うために必要な衛星画像の数です。

次に、作成した関数を使用して最新の 4 つの画像を取得します。


RY_latest, RY_latest_timestep = download_data()

画像を取得した後、wradlib の vis.plot_ppi メソッドを使用してデータをプロットします。

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

これでデータがロードされました。次にモデルをロードしましょう。

まず、関連するクラスをインポートします。この記事では TensorFlow を使用します。

tensorflow.keras.layers からインポート、Input、Conv2D、Activation、Concatenate、Conv2DTranspose、MaxPool2D

tensorflow.keras.models からモデルをインポート

3 つの基本的な構成要素を構築してみましょう。これらの「ビルディングブロック」は、この実装に従って、アーキテクチャ全体を作成するために使用されます。

最初のブロックは一連の畳み込み層に対応し、これを「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

2 番目のブロックは、エンコーダー部分 (畳み込みブロック + 最大プーリング) からの構築に使用されます。これを「encoder_block」と呼びます

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

3 番目の最後のブロックは「decoder_block」(アップコンボリューション + 連結 + コンボリューション) です。

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

これらの構成要素を組み合わせて 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

他の目的で使用したい場合は、自由に実装を試したり、改善したり、ニーズに合わせて調整したりしてください。

このモデルのトレーニングには多くの時間がかかる場合があります。幸いなことに、U-Net アーキテクチャに基づいて作成され、降水ナウキャスティングに特化した Rain-Net と呼ばれるモデルがあります。

GitHub リポジトリのクローンを作成しましょう

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

次に、このモデルの事前トレーニングされた重みをダウンロードします。

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

次のステップでは、リポジトリで見つかったアーキテクチャに基づいてモデルを作成し、このモデルにダウンロードされた重みをロードします。

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

ダウンロードした画像のサイズは 900*900 ピクセルです。 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

次に、予測を行う関数を作成します。

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

次に、前にダウンロードしたデータに対してこの関数を呼び出します

Y_pred = prediction(model, RY_latest)

予測をプロットし、保存した結果を使用して、予測を視覚化できる gif 画像を作成できます。

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

ここまで進んだことをおめでとうございます! Rain-Net を使用して予測を行い、それを視覚化できるようになりました。

＃＃結論

この記事では、機械学習モデル (Rain-Net) を使用して降水量ナウキャストを実行しました。 TensorFlow を使用して構築した U-Net アーキテクチャを使用しました。衛星画像を予測するために事前トレーニングされたモデルをロードしました。

この実装はさまざまな方法で改善できます。例えば：

データセット上のモデルを微調整する
アーキテクチャで CBAM (畳み込みブロックアテンションモジュール) などのアテンションベースのモジュールを使用します。

参考文献

無料ワークショップにぜひお越しください!

無料ワークショップでデータサイエンティストとしてのキャリアをスタートしましょう。ワークショップは適応可能なカリキュラムに基づいており、業界の専門家が指導します。