再解析データをプロットしてみよう！
寒冷渦ぽむぅ！坂浦いとです！
今回は，再解析データをプロットしてみようと思います．
虹ヶ咲学園 地学部気象班
上原歩夢，中須かすみ
環境
WSL Ubuntu 22.04.3
Python 3.10.12
Cartopy 0.22.0
matplotlib 3.5.2
numpy 1.22.3
pandas 1.4.2
xarray 2023.9.0
再解析データとは
かすみ：
歩夢先輩，再解析データとはなんでしょうか？
歩夢：
再解析データとは，過去の観測データを現在の気候数値モデルで同化したものだよ．
かすみ：
いわゆるデータ同化を用いているんですね
歩夢：
データ同化とは，シミュレーションの結果をたたき台にして観測値で補正することだよ．シミュレーションではコンピュータ内での数値計算であるため，現実の値と誤差が生じると考えられるよ．そのため現実の値である観測値で補正するよ．一方で観測値は観測点がある場所のみ定義されるためまばらに存在して，一方でコンピュータシミュレーションは空間内の格子点上で定義されていて，データ同化で観測と数値計算の長所を組み合わせて短所を埋め合わしていく感じだね．
かすみ：
再解析データはかすみんたちでもゲットできるんでしょうか？
歩夢：
再解析データは，ユーザー登録が必要なものとそうでないものがあるよ．今回はユーザー登録が不要であるNOAA(アメリカ海洋大気庁)が提供する「NCEP/DOE Reanalysis II」を使ってみよう．
かすみ：
ほかにはどういうものがあるんでしょうか
歩夢：
再解析データによって水平格子点間隔が違うけれど，今回使うNCEP/DOE Reanalysis IIは水平格子点間隔が2.5°×2.5°，ECMWF(ヨーロッパ中期予報センター)が提供するERA5は0.25°×0.25°，日本の気象庁が提供するJRA-55やJRA-3Qは1.25°×1.25°だよ．ほかにもあるけれど割愛するね.(単に作者が使ったことがないので無知)ファイル形式はNCEP/DOE Reanalysis IIやERA5はnetCDF形式(.nc)，JRA-55やJRA-3Qはgrib1形式だよ．今回は，ファイルサイズおよびPythonで読み取りやすい.nc形式であるNCEP/DOE Reanalysis IIをつかって気象データをプロットしてみよう．
データの入手とデータ構造
かすみ：
どうやって入手するか教えてください．
歩夢：
このサイトへ行ってみよう！
https://psl.noaa.gov/data/gridded/data.ncep.reanalysis2.html
ここではNCEP/DOE Reanalysis IIを配布しているよ．今回は1995年4月17日の天気図をプロットしてみよう．
そのため「hgt.1995.nc」をダウンロードして作業ディレクトリに置こう．
かすみ：
はい！このあとはどうすれば
歩夢：
Pythonを用いて読み込むよ．
じゃあこれを書いてみよう．
#!/bin/usr/env python3
###---使用モジュール---###
import xarray as xr # 多次元データを読むためのモジュールだよ
import numpy as np
import matplotlib.pyplot as plt # プロットに使うよ
import cartopy.crs as ccrs # 地理的データのプロットに必要だよ
import matplotlib.ticker as mticker # 緯度経度格子線のフォーマット調整に使うよ
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
# 緯度，経度のフォーマットの調整に使うよかすみ：
これをコピーするんですね！
歩夢：
走らせてみましょう！
かすみ：
やりましたよ？何もでないよ．
歩夢：
エラーが出ずうまくいったということだよ．
かすみ：
このあとどうすればいいでしょうか？
歩夢：
データを読んでみよう！
かすみ：
もう準備できたんですか？
歩夢：
.pyファイルと同じディレクトリにファイル置いているのを確認したら，これを書いてみよう．
(同じでもなくても大丈夫ですが，パスの設定が必要となります)
ds = xr.open_dataset("hgt.1995.nc") # xr.open_dataset()で.ncを読み込もう
print(ds) # 画面上に表示するよかすみ：
書きました！走らせます！わああなんか出てきました：
<xarray.Dataset>
Dimensions: (level: 17, lat: 73, lon: 144, time: 1460)
Coordinates:
* level (level) float32 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
* lat (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
* lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
* time (time) datetime64[ns] 1995-01-01 ... 1995-12-31T18:00:00
Data variables:
hgt (time, level, lat, lon) float32 ...
Attributes:
Conventions: CF-1.0
title: 4x daily NCEP/DOE Reanalysis 2
comments: Data is from \nNCEP/DOE AMIP-II Reanalysis (Reanalysis-2)...
platform: Model
source: NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Model
institution: National Centers for Environmental Prediction
dataset_title: NCEP-DOE AMIP-II Reanalysis
References: https://www.psl.noaa.gov/data/gridded/data.ncep.reanalysi...
source_url: http://www.cpc.ncep.noaa.gov/products/wesley/reanalysis2/
history: created 2002/03/15 by RHS (netCDF2.3)\nConverted to chunk...歩夢：
うまく読み取れたみたいだね！
かすみ：
でもこれって何言ってるんですか～
歩夢：
気候データについて学習しよう．
かすみ：
先に言ってくださいよ～
歩夢：
Datasetとはいろいろなデータ(Dataarray)の集合体と考えていいよ．
かすみ：
文字通りですね！
歩夢：
解読してみよう！まず，
<xarray.Dataset>
Dimensions: (level: 17, lat: 73, lon: 144, time: 1460)xarray.DatasetになったことでnetCDFがxarray.Datasetとして読み込まれたことを示しているね．
つぎにDimensionsとは，それぞれの次元に入っているデータの数だね．
かすみ：
「(level: 17, lat: 73, lon: 144, time: 1460)」とは，levelに17，緯度に73，経度に144，時間が1460あるということでしょうか？でもそもそもデータってどうなっているかイメージできません．
歩夢：
データはつぎのようになっているよ．
つまり，一般的に温度，降水量，高度などの変数がそれぞれ緯度，経度，高度，そして時間を引数として決まっているよ．数式で表すと，変数$${A}$$に対して，
$$
A = A(x,y,z,t)
$$
という関係を各時空点$${(x,y,z,t)}$$に引数数の次元でおいたのがDataarrayでDataarrayの集まりをDatasetと呼ぶよ．さらにcoordinates(座標系)は，軸$${(x,y,z,t)}$$ごとに定義されていて，たとえば$${x}$$では経度列が出てくるよ．
かすみ：
空間と時間を引数とする函数を座法軸上に収納した物なんですね．
歩夢：
そんな感じだよ．ここに関しては理屈よりかはたくさん使って慣れていこうね．ちなみに引数の順序は，$${(t,z,x,y)}$$となっていることに注意しようね．今回，気圧面だから$${z}$$は$${p}$$と書けるよ．
さてこちらを読んでみよう：
Coordinates:
* level (level) float32 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
* lat (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
* lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
* time (time) datetime64[ns] 1995-01-01 ... 1995-12-31T18:00:00これは座標値を示しているよ．
かすみ：
levelは気圧面の値に対応して，1000.0, 925.0,…,10.0と書いてありますね．
歩夢：
そうだよ．これが鉛直軸で，
かすみ：
latは緯度，lonは経度でたしかに最初に歩夢先輩が説明したとおり，2.5°刻みで全球で定義されていますね！
歩夢：
そうだね！そして時間は1995年1月1日から12月31日まで6時間ごとにあるよ！つぎに
Data variables:
hgt (time, level, lat, lon) float32 ...は，データ変数を示していて，hgtという変数は型がfloat32で，時間，気圧面，緯度，経度を軸とするよ！hgtとは高度であって単位は[m]だよ．さらに今回変数が一つだが，複数のことがあるよ．
かすみ：
あとは何でしょうか？
Attributes:
Conventions: CF-1.0
title: 4x daily NCEP/DOE Reanalysis 2
comments: Data is from \nNCEP/DOE AMIP-II Reanalysis (Reanalysis-2)...
platform: Model
source: NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Model
institution: National Centers for Environmental Prediction
dataset_title: NCEP-DOE AMIP-II Reanalysis
References: https://www.psl.noaa.gov/data/gridded/data.ncep.reanalysi...
source_url: http://www.cpc.ncep.noaa.gov/products/wesley/reanalysis2/
history: created 2002/03/15 by RHS (netCDF2.3)\nConverted to chunk...歩夢：
これはAttributes，つまり属性で今回つかうデータの情報が書いてあるよ．でもあまりデータ解析にはあまり関係ないよ！
かすみ：
そうなんですね！
高度場の入手
歩夢：
さて今回は高度データを使いたいから，こう書いてみよう：
hgt = ds.hgt # 「xarray.Dataset.変数名」 で変数のdataarrayを抽出できるよ．こうすることで，Datasetに入っている変数を抽出できるよ．
変数hgtを抽出してhgtに収納しよう．
一旦どういうデータか表示するよ：
print(hgt)<xarray.DataArray 'hgt' (time: 1460, level: 17, lat: 73, lon: 144)>
[260907840 values with dtype=float32]
Coordinates:
* level (level) float32 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
* lat (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
* lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
* time (time) datetime64[ns] 1995-01-01 ... 1995-12-31T18:00:00
Attributes: (12/13)
long_name: 6-hourly Geopotential Heights on Pressure Levels
units: m
precision: 0
GRIB_id: 7
GRIB_name: HGT
var_desc: Geopotential height
... ...
level_desc: Pressure Levels
statistic: Individual Obs
parent_stat: Other
standard_name: geopotential_height
actual_range: [ -572. 32265.]
valid_range: [-1500. 35800.]かすみ：
またなんか出てきました！
歩夢：
解読してみよう．
<xarray.DataArray 'hgt' (time: 1460, level: 17, lat: 73, lon: 144)>
[260907840 values with dtype=float32]これは，変数hgt(高度)について，軸は時間，気圧面，緯度，経度になっていることだよ．値の数は260907840個でデータ型はfloat32だよ．
かすみ：
多すぎます～！
歩夢：
つぎに
Coordinates:
* level (level) float32 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
* lat (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
* lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
* time (time) datetime64[ns] 1995-01-01 ... 1995-12-31T18:00:00は，座標軸だよ．さっきと同じだね．
かすみ：
あとの部分は属性ですよね？
歩夢：
そうだよ．さて，今回は1995年4月17日0:00(UTC)の500hPaのデータを図示したいから，これを入れるよ：
hgt = hgt.sel(time ="1995-04-17T00:00", level = 500)
# xarray.Dataarray.sel(取り出すcoordinateの値，複数指定OK)
print(hgt)これによって，1995年4月17日0:00の500hPaのデータが得られたよ．走らせると，
<xarray.DataArray 'hgt' (lat: 73, lon: 144)>
[10512 values with dtype=float32]
Coordinates:
level float32 500.0
* lat (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
* lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
time datetime64[ns] 1995-04-17
Attributes: (12/13)
long_name: 6-hourly Geopotential Heights on Pressure Levels
units: m
precision: 0
GRIB_id: 7
GRIB_name: HGT
var_desc: Geopotential height
... ...
level_desc: Pressure Levels
statistic: Individual Obs
parent_stat: Other
standard_name: geopotential_height
actual_range: [ -572. 32265.]
valid_range: [-1500. 35800.]かすみ：
levelが500になって，timeが1995-04-17になりましたね！
歩夢：
もうひとつ重要なことがあって，さっき次元が(時間，気圧面，緯度，経度)だったけれど次元が削減されて(緯度，経度)の二次元になったよ．これによって天気図がプロットできるね！
かすみ：
ところで天気図ってどうやってプロットするんですか？
気象データのプロット
歩夢：
天気図をプロットするためにはmatplotlib.pyplotという描画ツールと地理的データを表現するのに必要なcartopyなどを用いるよ．描画は設定が難しいから，これを写経してね：
def plot_contour(x, y, data, extent, limits, annotation, interval, title):
fig = plt.figure(figsize=[16, 9]) # figオブジェクトを用意(サイズは横16cm×9cm)するよ
ax = fig.add_subplot(projection=ccrs.PlateCarree(int((extent[0]+extent[1])/2)))
# axオブジェクトを用意，メルカトル図法でプロットするよ
ax.set_extent(extent, ccrs.PlateCarree())
# 描画範囲extentに絞るよ
# plot basemap
ax.coastlines(resolution="10m") # 海岸線を描くよ．
gl = ax.gridlines(crs=ccrs.PlateCarree(), linestyle=":", draw_labels=True)
# 格子線を用意するよ．線のスタイルは点線でラベルも書くよ
gl.top_labels = False # 北縁のラベルは書かないよ
gl.right_labels = False # 東縁のラベルは書かないよ
gl.xlocator = mticker.FixedLocator(np.arange(-180,180,20)) # 20°間隔
gl.ylocator = mticker.FixedLocator(np.arange(-90,91,10)) # 10°間隔
# 緯度経度ラベルを書くところを決めるよ．
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())
# 緯度，経度軸の記法を定まるよ
anno = np.arange(limits[0], limits[1], annotation)
# 値付きコンタの系列を決めるよ
inter = np.arange(limits[0], limits[1], interval)
# 値なしのコンタの系列を決めるよ．
ax.contour(x, y, data, # コンタ図のプロット: 経度データ(1D)，緯度データ(1D)，二次元データ
levels=inter, # コンタが引かれる場所
linewidths= 0.3, # 線の太さ
colors = "black", # 線の色
transform=ccrs.PlateCarree()) # 地理データをプロットするのに必要
clabels = ax.contour(x, y, data,
levels=anno,
linewidths= 0.5,
colors = "black",
transform=ccrs.PlateCarree())
# 値付きコンタのプロット，オブジェクトをclabelに収納
ax.clabel(clabels) # 値を入れる．
ax.set_title(title) # タイトルを入力するよ
plt.show() # 画面上に描画この函数をプロットしたら早速図示してみよう：
かすみ：
でも，x軸，y軸はどうやって出すんですか？
歩夢：
こうやって取得するよ
x = hgt.lon
y = hgt.lat
# 変数.座標軸で座標軸のdataarrayを入手できるよかすみ：
座標軸もデータ扱いするんですね．わざわざ自分でnp.arangeやnp.linspaceで座標軸を作らなくてよいんですね！
歩夢：
そこがxarrayないしnetCDFの強みだよ！
さてプロットしてみよう．
plot_contour(経度，緯度，2Dデータ，範囲[最小経度，最大経度，最小緯度，最大緯度]，変数の値域，値付きコンタ間隔，値なしコンタ間隔，タイトル)
で
plot_contour(hgt.lon, hgt.lat, hgt,
[0, 360, -90, 90], [0, 8000], 120, 40,
"500hPa 00:00 17 Apr 1995") かすみちゃん，押してみて！
かすみ：
わかりました！出てきました：
歩夢：
これが1995年4月17日の全球の500hPa高度場でよ．中緯度では，等高線が波打っていて偏西風波動に対応するよ．一方，熱帯では等高線が少ないね．
かすみ：
日本付近に注目してみたいです！
歩夢
plot_contour()のextent(描画範囲)を次のように変えてみよう！
plot_contour(hgt.lon, hgt.lat, hgt,
[100, 180, 15, 70], [0, 8000], 60, 20,
"500hPa 00:00 17 Apr 1995") ちょっと，コンタ間隔変えてみたよ．
かすみ：
結果出しますね！
出てきました！
楽しいですね！
歩夢：
4月17日の日本付近では，西日本ではゾーナルで，オホーツク海にはブロッキングリッジがみられますね！アリューシャン列島には寒冷渦ぽむぅかな？
かすみ：
いままで，気象庁など探していろいろデータ探していてもなかなか過去の天気図って手に入れられにくいし任意の気圧面をプロットしたいと思っていたけれどついに自由にプロットして天気図検討できるようになりましたね！
歩夢：
再解析データを使うことによって大きな気象イベントがあった日やかすみちゃんの誕生日の天気図を作り出すことができるよ！今回は気圧面の高度面をプロットしたけれど，温度や風，海面更正気圧や地上気温もできるよ！さらにPythonを使っているから相関解析や回帰分析，EOF(経験的直交函数)分析，いわゆる主成分分析もできて自由度がすごく高いよ！
かすみ：
可能性がかなり広がりましたね！今度教えてください！
歩夢：
今度教えるね！おつかれさまかすみちゃん！
かすみ：
歩夢先輩，ありがとうございました！
問題
かすみんレベル：
別の気圧面，別の時間のデータをプロットしてみよう．
とくにあなたの誕生日の天気図を作ってみよう．
ぽむぅレベル：
気圧面データで別の変数のデータを入手してプロットしてみよう．
りなりーレベル：
地上天気図を作成してみよう．
参考文献
気象データ解析 - Tsubasa Kohyama https://sites.google.com/site/tsubasakohyama/%E6%B0%97%E8%B1%A1%E3%83%87%E3%83%BC%E3%82%BF%E8%A7%A3%E6%9E%90
気象データ解析のためのmatplotlibの使い方 (作成者：山下陽介（国立環境研究所)) https://yyousuke.github.io/matplotlib/index.html
xarray公式ドキュメント https://docs.xarray.dev/en/stable/