氷河合宿その５f：気圧と標高の関係

 月平均データについては、23層200 hPaまでのダウンロードを推奨しています。

が、ヒマラヤみたいに9000 m 近くまで地面があるところならともかく、グリーンランドなら3000 mちょいしかないわけで、「そんなに高いところまでいらないんじゃね？」と思うのは当然です。層を減らせばダウンロードにかかる時間もデータ容量も減らせますしね。

そこで、配布しているパッケージに含まれる、「z_pressure.py」の出番です。

自分のサイト用に必要なところをいじって走らせると、「zp_【サイト名】.csv」というファイルが生成されます。これは、対象期間中のターゲットサイトを含むグリッドにおける、気圧と最低高度、平均高度、最高高度、を示しています。

グラフにしてみるとこんな感じ。横軸が気圧、縦軸が高度です。自分が扱うサイトのおおよその標高を把握しておけば、何 hPa までをダウンロードすればよいか、当たりをつけられると思います。hourly データのダウンロードを始める前に確認しておくとよいでしょう。例えば、グリーンランドは4000 m までカバーしておけば十分でしょうから、550 hPa まででよい、とか。
※hourly の等気圧面ジオポテンシャル高度は月平均よりバタつくので、余裕をもってダウンロードしておくことをお勧めします（気圧低め、高度高め）。

スクリプトの中で工夫したのは、numpy.mean / numpy.min / numpy.max で各年の平均、最小、最大を求め、出力するときにさらに平均、最小、最大を求めているところです。

氷河合宿その５e：ポイントデータの抽出（等気圧面月平均）

地表面データだけでも観測データとの比較や質量収支モデルを走らせることは可能ですが、今回紹介する等気圧面データからの気温抽出は、ウチのオリジナルアイデアで、かつ、今のところ他の気温推定方法よりも成績が良いという実績があります。

今回のパッケージはこちらからダウンロードしてください。
【補足】「危ないサイトにはアクセスできません」的なメッセージが出た場合、「上のアドレスに表示される「https://」から s を削除するとダウンロードができるみたいです。（https へ移行しなきゃ、と思いつつサボってますすみません）

Pythonで「pl_extract_month.py」を実行すると、等気圧面データから抽出した気温（pt）と気温減率（lap）、その抽出に使ったデータが生成されます。

やっていることは、ターゲットのサイトの標高を挟む等気圧面をジオポテンシャル高度から探しだし、その二つの高度と気温から気温減率を求め、さらにターゲット標高の気温を求めています。

tl、zl、pl：ターゲットより低い側の等気圧面における、気温（tl）、高度（zl）、気圧（pl）

tu、zu、pu：ターゲットより高い側の気温（tu）、高度（zu）、気圧（pu）

例えば、サンプルデータを実行してできるSigma-Bサイトの高度（zl、zu）をグラフにしてみると、

ターゲット標高の944 mを挟んで変動していることがわかります。また、その気圧面（pl、pu）はどうかというと、

この高さ周辺での気圧面は925、900、875で、それらが変化していることがわかります（pl の方が高い値）。

対面合宿時には、前回抽出した地表面データのオリジナル気温 t2m とそれを気圧減率 -6.5 ℃/km で補正した気温 at と一緒に観測データと比較した結果を紹介してもらいます。

スクリプトの解説

前半の大部分は地表面データに対する処理と同じです。

インポートしてサイト情報を読み込み

ヘッダーを書き込んでおいて、年ごとに気圧面データを読み込み、変数を設定。このデータは高さ方向のデータ（nz と dpl）も含んでいるのが地表面データと違うところ。

ターゲットサイトの座標を決めるやり方とそのコツは前回解説しました。

そしてこれ以下の部分で「ターゲットサイトの標高を挟む二つの高度」を探しています。izを 1 から一つずつ増やしていって（L60）、もし、上回ったら（L53）、上側（izu）と下側（izl）を決めてやり、izl の拘束（= nz）が外れることで while ループが終了するようになっています。

鋭い人は、「Python はゼロ始まりじゃなかったっけ？」と思うかもしれません。ここでは、「上側の等気圧面がターゲット標高を上回っているか？」を探しているので、1 から始めています。

もし、0 から始めていきなりターゲット標高を上回っていたら、izl が ‐1 になってしまって計算ができなくなりますよね？※こういうケースでは、1000 hPa と 975 hPa から外挿して求めていることになります。

二つの等気圧面が決まったら、次の時間ステップのために iz を二面ほど下げておきます（L58）。次の月に気圧が下がる（ジオポテンシャル高度が上がる）と、ターゲット標高がいきなり低い側の等気圧面高度よりも低くなり、内挿でなく、外挿で求めてしまう可能性があるので、常に「より低いところ」から探し出すようにするためです。※でも、ヒマラヤとか5000 mもあるようなところで、毎回1000 hPaから探し出すのは非効率

L56 では、この下げ幅が十分でなかった場合にアラートが出るようにしています。

L61-65 では、ダウンロードの設定ミスで、ターゲット標高がデータセットの最高高度よりも高かった場合、とりあえず外挿で値を決めつつ、アラートを出しています。

あとは、izu、izl の気温、標高から気温減率を求め（L72）、ターゲット標高の気温を決めています（L73）。二次元上の二点間の傾き（＝気温減率）と、その間にある点の値を決めています。

氷河合宿その５d：抽出グリッドの決め方

【訂正】グリッドサイズを一桁間違えてた

抽出したいサイトを含むグリッド（の 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＝一番目」にあると認識されなければいけません。という端っこの処理をしているのがこちらです。

「もし東西方向のデータ数以上になったら、減らしてね」という感じ。この境界に氷河があるかどうかも知らないけど、まぁちょっとしたこだわりポイントでした。

氷河合宿その５c：ポイントデータの抽出（地表面月平均）

ダウンロードがひと段落したら、いよいよ自分が解析したいサイトの情報を取り出すことに挑戦します。

ダウンロードがうまくいってない人のために、北西グリーンランドあたりのデータを用意しました。こちらからダウンロードしてください（クリックしたらダウンロードされるはず）。
【補足】「危ないサイトにはアクセスできません」的なメッセージが出た場合、「上のアドレスに表示される「https://」から s を削除するとダウンロードができるみたいです。（https へ移行しなきゃ、と思いつつサボってますすみません）

これを作業フォルダで展開してください。※単にクリックするだけだと中身を閲覧してるだけで、スクリプトは動きませんからね。

フォルダ内でターミナルを立ち上げ、ipython からの run met_share_month.py ですぐに動くはずです。（具体的なやり方は知りませんが）jupyter や spyder でも動くはず（たぶん）。まずは動くかどうか試してみてください。

自分のサイト用にするには

まず、site.ini をテキストエディタで開きます。自分でダウンロードが済んでいる人は、欲しいところの座標と標高情報、期間を指定します。よくある間違えは経度と緯度を取り違えることです。西経はマイナス表記にしてください。南半球は緯度をマイナス表記にします。標高は実数表記にしてますが、ここは整数でも大丈夫。元のnetcdfファイルの範囲に含まれているなら、複数のサイトを一括して抽出することが可能なようにしてあります。もちろん一カ所でもOK。一行目のヘッダーは消さないこと。

met_share_month.py の L8-9 のフォルダとファイル名のルールを自分のモノに合わせる

だけでいけるはず。

スクリプトの説明

以下、met_share_month.py で何をやっているかを解説します。

最初の数行はパッケージの導入です。これでエラーが出るようなら、パッケージをインストールしてください。netCDF4以外はpaythonと一緒に入っていると思いますが。

L8-9　データが入っているフォルダとその名前を指定

L11-L15　site.iniからターゲットとするサイトの個数（L12）を読み取り、座標情報をそれぞれの配列に移し替える。

L16-L21　出力ファイル名を設定。サイトの数だけ用意し、それぞれのファイルにヘッダーを書き込む。（この「空のテキスト配列を設定する」やり方見つけるの苦労した）

L23-L33　ERA5の地表面標高を読み込む。dlon, dlat は経度緯度情報をマルっと代入。delv は標高の配列で、抽出サイト分だけ用意。nx, ny は東西と南北のデータ数、dgx, dgy は東西方向、南北方向それぞれの空間解像度、lx は東西方向の西の境界（x軸の左端、という意味合い）、uy は南北方向の北の境界（y軸の上端、という意味合い）

L34-36　まとめファイルを開き、ヘッダーを書き込む。

L37-　サイト毎に対応グリッドを探す。

L38-41　以前説明したとおり、全球をカバーする場合は経度が0-360表記なので、site.ini の経度が西経マイナス表記の場合、360を足している（L37）。L38でgeopotentialのグリッドで左から何個目に相当するかを計算。

L42　緯度方向に北から何個目にあるかを計算

L43　標高に換算した上で配列に代入

L44-46　確認のため画面にサイトとgeopotential の緯度経度標高を表示しつつ、ファイルに書き込む

L48-49　対象期間を指定、ダウンロード済みのデータがあることが前提

L49　暫定的な気温減率を設定

L51-63　東西南北のデータ数を確認したり、解像度、境界を設定（上のL23-33と同じ）

L65-　サイトごとに処理

L66　site.ini で指定した解析対象期間に該当するかどうかの判定

L67-70　東西南北何個目のグリッドに該当するかをチェックし、NetCDFデータの範囲外なら計算を中止する

L72-　ファイルへ出力、それぞれのパラメータについて、単位を変換したりしてる

L80-87　「相対湿度」というパラメータが無いので、露点温度と気温から求めている

L89-93　降水量は一日あたりの値なので、月降水量に変換するには各月の日数が必要なので、それを計算

L96-97　ファイルに書き込み

という具合です。

この中で、西と北の境界（lx, uy）とグリッド何個目かを計算しているところは工夫のしどころなので、別に解説します。

氷河合宿その４h：ダウンロードが進まない！

hourlyデータのダウンロード、合宿参加予定者から「一晩たっても全然ダウンロードできないのですが、こういうもんですか？」との悲鳴が来ています。かくいう自分も全然ダウンロードが進まず、「合宿までに終わるのか？」と不安になっていました。

今朝方、四月から学振PDとしてきているKKさんが有効そうな解決策を発見、教えてくれました。

ERA5データのあるCDSにいき、ログインしてから「Your Request」をクリック

待機中やダウンロード済みのもの全て選択して片っ端から全削除

でもって、も一度ダウンロードを再開します。

この時、スクリプトを書き換えて、既にあるファイルを再度ダウンロードしないようにしましょう。

追記

KKさん調べで、スクリプト内の

client = cdsapi.Client()

を

client = cdsapi.Client(delete = True)

としておくと、ダウンロード毎にERA5サイト内をお掃除してくれるので、待機時間が長くならない模様

配布用のスクリプトに反映しておきました（リンク先の記事の一番下）。

氷河合宿その４g：ダウンロードスクリプト応用編

 合宿参加者が提出してくれたスクリプトの中で、「おっ、これは使える」と思ったモノを紹介します。

こちらはフォルダの有無を確認し、なければ作成してからダウンロードを開始するというもの。

import os

odir = 'nc'

os.makedirs(odir, exist_ok=True)

        |

target = f'{odir}/site{yr}{mn}mon.nc'

こちらは変更しがちなパラメータについて、冒頭で指定しておくというもの。

syr, eyr = 2010, 2020

odir = 'nc'

site = 'arctic'

latn, lats, lonw, lone = 80, 70, 60, 100 # 北南西東で指定

        |

"area": [latn, lonw, lats, lone] # ERA5の指定に合わせて「北西南東」となっているのに注意

target = f'{odir}/{site}{yr}{mo}hour.nc'

#もしくは 

target = odir + '/' + site + yr + mn + 'hour.nc' 

それぞれ工夫してくれたおかげで自分も新しいことも覚えられてなかなか良かったです。

氷河合宿その４f：ERA5 hourlyデータ連続ダウンロード（他、ダウンロード用スクリプト）

前回出した課題の答えです。

こちらからお手本スクリプトを入手してください（右クリックから「名前を付けて保存」で入手）。

以下、内容の説明です。

for iy in range(2012,2014,1):

の中に、

for im in range(1,13,1):

として、月ループを用意します。12ヶ月ある場合、13まで指定しないといけないのはpythonのお作法。これ以降、最後までインデントをもう一つ下げる必要があります。

で、次の行がポイントで、

mn = format(im,'02')

ここで、数字が入っている im を、文字列の mn とし、かつ、二桁指定で不足する場合はゼロをつけています。

「python ゼロ埋め」

で検索すると、やり方はいろいろあるようで、

mn = str(im).zfill(2)

mn = str(im).rjust(2, '0')

mn = f'{im:02}' # mn = f'{im:02d}' もあり

mn = '{0:02}'.format(im)

mn = '%02d' % im  # 古いやり方のようで推奨されてないみたい

でも、同様の処理ができるようです。Pythonのバージョンによっては動かないのもあるらしい。

試しに im = 2 に続けて上のやり方をそれぞれ試し、im リターン、mn リターンとして中身を表示してみましょう。

imが

2

mnが

'02'

と表示されれば、mnが二桁表記の文字列であることを示しています。

さらにベタなやり方は

 

if im < 10: # 一桁かどうかを判定

    mn = '0' + str(im) # 一桁の場合、前にゼロを追加

else:

    mn = str(im) # 二桁ならそのままテキストに

 

というものです。他にも、

 

mn = str(im) # まずテキストにする

if len(mn) == 1: # len()は配列の長さをみる

    mn = '0' + mn # 一桁ならゼロ追加

というのもありました。なかなか工夫していますね。

 ただ、二桁ならまだいいですが、三桁、四桁と増えていくとこのやり方ではちょっと面倒ですね。

以降はmonthlyと似たような構造です。

"year"でyr、"month"でmnが指定され、

"day" で31日分、"time"で24時間分が指定されていることがわかりますね？

あとはファイル名のところで月を追加しているのが変更点です。

dayやtimeをループを作って書いていた猛者もいましたが、それがうまく機能するかは確証がありません。ERA5サイトで作成されるAPIスクリプトをコピーしてきて、yearとmonthのところだけを置き換える方が安全です（インデント下げるのを忘れないように）。

以下に、自由大気（pressure level）のmonthly, hourlyのデータダウンロードのスクリプトも併せて置いときますので、各自好きなように使ってください。氷河合宿に関係のない人もご自由に

pressure levelのhourlyデータはの高さ方向については、自分が必要なサイトの標高のおおよその気圧がわかっている人は、-150hPaくらいまでを入手しておけばよいです。提供スクリプトの場合23層あるので、できるだけ狭い範囲を指定することをお勧めします。スンゲー時間かかる。

なお、monthlyとhourlyで範囲やlevel数をそろえておく必要はありません。srfとplの範囲もバラバラで大丈夫（ターゲットのサイトが含まれている必要はある）。

get_srf5mon.py # 地表面monthlyデータ

get_srf5hour.py # 地表面hourlyデータ

get_pl5mon.py # 自由大気monthlyデータ（気温とジオポテンシャル高度）

get_pl5hour.py # 自由大気hourlyデータ（気温とジオポテンシャル高度）

get_srf5geopotential.py # 地表面ジオポテンシャル高度