【訂正】グリッドサイズを一桁間違えてた
抽出したいサイトを含むグリッド(の x と y )をどう決めるかについて解説します。
ERA5のデータは、この図の左側のような配置になっています。
例えば、経度0°のグリッドには、-0.125°~+0.125°が含まれています。
ちなみに、NCEPは右のような配置
一つのやり方としては、各グリッドの中心座標とターゲットサイトの座標の距離を計算し、その距離が最も小さいグリッドを選ぶ、というやり方があります。ただ、全球データを扱う場合、1440 + 721 グリッドを確認しなければならず(dx が一番小さい and dy が一番小さい)無駄が多いです(まぁ最近のPCは速いのですぐ終わるんですけどね)。
ここで利用するのが、int( ) です。これは、実数を整数に変換する関数です。0.1 ずつ増える実数をint( ) で整数にするスクリプトを試しに動かしてみましょう。
0~0.9 は 0、1.0~1.9 は 1 になっていることがわかります。これが曲者なのは、マイナス方向にも同じように働くので、-0.9~0.9 が「0」になることです。なので、ゼロをまたぐ領域では注意が必要です。
スクリプトの中では、「西の端(データの左端)から何個目か」「北の端(データの上端)から何個目か」を算出しています。
で、この「端」を設定する際、グリッドサイズの半個分、外側にずらしているのがポイントです(西端は西側へマイナス、北端は北側へプラス)。こうすることで、±0.125° を0~0.25°という範囲になるようにしています。
あとはグリッドサイズ(0.25°)で割ってやって、整数にする関数をかませば、グリッドで何個目か、という値が得られます。
なお、緯度方向は「北から南」に並んでいるので、緯度変数にマイナスがかかっています。
自分で用意する際は、プラスやマイナスが正しく設定できているか、確認のために求まった緯度経度とサイトの緯度経度を画面に表示させると良いでしょう。両者の差が0.125° 以内に収まっていれば合っている、ということになります。
以降は興味のある人向け:超細けえ話
ところで、全球データの場合、経度は0~360°です。このデータと突き合わせる場合、西経のマイナスは360を足して180~360になるように変換しています。この時、0°ギリギリの359.875~360°にあるサイトはどこに含まれるでしょうか?
計算してみるとわかりますが、(西端の境界が0.125西側にずれているので)1440という値が得られます。「ERA5のグリッド数は360 / 2.5 = 1440だからOKじゃね?」って思いますよね?ところがPythonでは、1440の配列は、0~1439です。だから、このままこのグリッドを呼び出そうとすると「ねーよ」とエラーがでます。
この部分は本来は「0=一番目」にあると認識されなければいけません。という端っこの処理をしているのがこちらです。
No comments:
Post a Comment