Monday, 18 November 2024

氷河合宿その5f:気圧と標高の関係

 月平均データについては、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 で各年の平均、最小、最大を求め、出力するときにさらに平均、最小、最大を求めているところです。

氷河合宿その5e:ポイントデータの抽出(等気圧面月平均)

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

今回のパッケージはこちらからダウンロードしてください。
【補足】「危ないサイトにはアクセスできません」的なメッセージが出た場合、「上のアドレスに表示される「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)。二次元上の二点間の傾き(=気温減率)と、その間にある点の値を決めています。


Friday, 15 November 2024

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

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

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


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

氷河合宿その5c:ポイントデータの抽出(地表面月平均)

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

ダウンロードがうまくいってない人のために、北西グリーンランドあたりのデータを用意しました。こちらからダウンロードしてください(クリックしたらダウンロードされるはず)。
【補足】「危ないサイトにはアクセスできません」的なメッセージが出た場合、「上のアドレスに表示される「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)とグリッド何個目かを計算しているところは工夫のしどころなので、別に解説します。

Thursday, 14 November 2024

氷河合宿その4h:ダウンロードが進まない!

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

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

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

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

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

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


追記

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

client = cdsapi.Client()

client = cdsapi.Client(delete = True)

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

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

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

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

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

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' 

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

氷河合宿その4f: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 # 地表面ジオポテンシャル高度

Tuesday, 12 November 2024

氷河合宿その4e:ERA5 pressure levelデータ

ダウンロードに掛かる時間などを考えると、来週の対面合宿に間に合うか心配になってきたので、課題とは別にダウンロードしてもらうことにしました。

今回は大気の等気圧面における気温とジオポテンシャル高度のデータです。ジオポテンシャル高度とは、指定された気圧面の標高のことです。
※入っている数値自体は標高ではなく、ちとややこしい。

地表面データだけでも質量収支モデルは動かせますが、ERA5の中で使われている標高と、サイトの標高にはずれがあります。このずれを補正するには気温減率が必要になりますが、地表面データには含まれていません。

また、ここがウチのオリジナルアイデアなのですが、等気圧面の気温と標高から、サイト標高を挟む二つの等圧面を探し出し、その間を線形近似することで、サイト標高の気温とその時その標高付近の気温減率を求める、という処理を行います。対面合宿&このブログで今後紹介しますが、これまでのいくつかの研究で、この方法で推定した気温が観測された気温と最もよく合うことがわかっています。

というわけで、こちらからスクリプトを入手してください(右クリック→名前を付けてリンク先保存)。以下説明


もう年ループは大丈夫ですね?

データセットやプロダクトタイプは指定し忘れないように。今回は気温とジオポテンシャル高度ですが、他にもいろいろパラメータがあります。


自由大気の場合、気圧面=pressure_levelを指定してやる必要があります。地面の一番高いところは8850mで、250~300hPaなのと、この後気圧と標高の関係を見るのに使うので、200hPaまで入手しておいてもらえると良いと思います。


後は地表面の月平均データと同じ。名前が重複してると上書きしてしまうので要注意。


氷河合宿その4d:ERA5月平均データ連続ダウンロード

単発のデータダウンロードができたら、次は複数年の自動ダウンロードをやります。

研究室サイトにスクリプトを置いたので、右クリックから「名前を付けて保存」で入手し、テキストエディタで開いてください。(若干修正)

三行目の

for iyr in range(2010,2021,1):

は、Pythonの「繰り返し」コマンドで、iyrというパラメータに、2010から2020まで、1ずつ増やして、という意味です。ん?2021までじゃないの?って思いますよね、Pythonはn個の行列があった場合、0から始まってn-1で終わる、というルールになってます。なので、2021と書いてある場合、2020まで、ということになります。

で、さらに重要なのが、このfor文で繰り返すコマンドは、以降の段落(インデント)が一つ下がっているコマンド全てになります。試しに、テキストエディタ上でfor文最後の:の後ろでリターンキーを押してみると、自動的にインデントが一つ下がったところにカーソルがあることがわかります。

yr=str(iyr)

print(yr)

 iyrは整数なので、次のstr(iyr)でテキストとして割り当て、画面上にその年を表示させています。(いつのデータを待機中なのかがわかるように)

以降しばらくはお試しダウンロードで作ったスクリプトと同じはずです。(ただし、年の指定が置き換わっているのと、インデントが下がってる)

ちゃんと動くスクリプトが用意できたら(今回のはまさにそうなのですが)、最後の"area"と保存ファイル名を変更するだけで、他の領域のダウンロードが可能になります。(欲しいパラメータが違ってくるなら作り直す必要ありですが)

ちなみに領域指定は

"area": [80, -80, 70, -60]

"area": [北, 西, 南, 東] ですからね。

最後のtargetで、保存ファイル名を指定していますが、

target = 'nc/arctic' + yr + 'mon.nc'

「nc/」が保存先のフォルダで、その中にarcticYYYYmon.ncというファイル(YYYYが年)が保存されるようにしています。

あとは、作業フォルダの中にncフォルダを作成した上で、ipythonを起動し、

run get_srf5mon.py

を実行すればダウンロードできるはずです。お試しあれ。

なお、全球をダウンロードすると、一年あたり約250MBになります(10年で2.5GB、40年で10GB)。許容範囲の大きさですが、先の記事に書いた通り、範囲指定と全球では、西経の扱いが違うので、以降の演習で提供する抽出スクリプトが使えないのと、ダウンロードの待機に時間がかかるので、推奨しません。

なお、Sublime Textでは、インデントを変更したい行を選択し、「Edit(編集)」→「Line(行の操作)」→「Reindent(インデント最適化)」で一括でインデントを下げることができるので便利。

合宿参加者への課題:hourly dataのダウンロードスクリプト作成

質量収支モデルに使用する、hourly dataを連続ダウンロードするスクリプトを作成しましょう。monthly data(範囲×要素×12ヶ月)に対して、hourly dataは、「範囲×要素×24時間×365日」と一気に730倍のデータ量になります。大量のデータリクエストでCDSサイトの反応がとても遅くなることもあり、月毎にダウンロードすることを推奨しています。(また、この後使う、ポイント抽出スクリプトは月毎にファイルが作成されていることを前提に作ってあります)

というわけで、月平均データのための上記のスクリプトに対し、datasetとrequestの中身をhourly dataに合わせることに加え、「月ループ」を作る必要があります。

注意点
参照サイトが違うので、APIスクリプトを用意してみて、どこが違っているか確認しておきましょう。
年ループの中に、更に月ループを用意すること
月を指定するとき、1, 2,,, 11, 12ではなく、01, 02,,, 11, 12とする必要があること

この二点を自分で調べて作成してみましょう。2日後(11/14木)にお手本スクリプトを解説付きで公開します。

なお、月の中の日数指定は、「select all」で31日分全選択すればいいです。31日に満たない月は、データがある分しかダウンロードされない。

Monday, 11 November 2024

氷河合宿その4c:ダウンロードトラブルあれこれ

 先週末から受講生から届いた「ダウンロードでけへん」問題をここに掲載しますので、これから取り組む人は参考にしてください。

※プロダクトの指定忘れ

逆に驚いたのですが、プロダクトの指定を忘れてもダウンロード自体はできるみたいです。ただ、pythonで中身を確認しようとしても、ファイルを読み込む時点でエラーが出ます。


※スクリプトの誤り

ERA5のサイトで生成されるスクリプトの、最後の「.download」の消し忘れでダウンロードがうまくいかないエラーが発生します。(ダウンロード直前まではちゃんと動くのに)

ちなみに、すぐに対応して再度トライすると、生成されたncファイルがERA5のサーバーに残っているので、(ちゃんと直せていたら)すぐにダウンロードが始まるはずです。


基本的には「4a:ダウンロードしてみよう」をちゃんと読みなはれ!という話なので、エラーが出たら、エラーメッセージをよく読むことと、元のチュートリアルを忠実になぞっているかどうか、今一度確認してみましょう。

追記
※「.cdsapirc」の間違い

氷河合宿その3a:ECMWF提供のコードをPCに登録」で、

「.cdsapirc」というファイルに

url: https://cds.climate.copernicus.eu/api
key: XXXXXXXX-XXXX-XXXX-XXXX-XXXXXXXXXXXX

をコピペして保存します。

ということを書きましたが、この際に、url: と key: を入れずにコピペしたことで、ダウンロードできない人が複数名いました。なぜわざわざ消す!!


Thursday, 7 November 2024

氷河合宿その3c:ECMWF提供のコード(.cdsapirc)の登録(トラブル解決編)

 お試しダウンロードがうまくいかない事例があり、解決したので対処法を紹介しておきます。

PCの構成:Windows11にWSL/Ubuntuをインストール

関連するブログの記述を参考に、home/user/に「.cdsapirc」を置く。

で、home/user/DLDLというフォルダにダウンロード用のスクリプト「test.py」を置いて

ipython

run test.py

を実行したものの、次のエラーが出てします。


途中、client = cdsapi.Client() でエラーが生じていることがわかりますが、大事なのは一番最後のメッセージで、

Missing/incomplete configuration file: /root/.cdsapirc

とあります。つまり、rootフォルダの中に「.cdsapirc」が無いよ、というメッセージ。home/user/に置いてあった「.cdsapirc」をrootにコピーし(警告が出たけどスルー)、再度Ubuntuから立ち上げ直して実行したら、上手くダウンロードできました。

エラーが出たら、メッセージをちゃんと読め、ということですね。


Tuesday, 5 November 2024

氷河合宿その5b:経度(と緯度)のルール

 引き続き、NetCDFの中身を確認していきます。

nc['longitude'][:]

と打ってみると、経度情報がガーっと表示されます。[:] は、「データ全部」という意味合い。


では、グリーンランド用にダウンロードしたファイルに対して同じようにしてみると。

西経がマイナス表記になっていることがわかりますね。


さらに、標高データとして全球分をダウンロードしたgeopotential.ncについても同様に確認してみましょう。


何かおかしいことに気づきましたか?

そう、西経がプラス表記で+180~+360で表現されています。(なお、データ数が多いと途中が省略されるみたい)

部分ダウンロードと全球ダウンロードでは、開始経度が違うことと、西経の表記が違うのです。


これの何が問題か?

この後のデータ処理では、データを抽出したい地点の緯度経度が、ERA5データのどのグリッドに含まれているか、を計算します。全域をソートして最も距離が近い点を選択、というやり方もアリなのですが、コードが長くなるのと、多点を抽出したい場合や、全球データの場合、時間がかかるのが難点です。

ともあれ、この後紹介するPythonのコードでは、それらに対応した処理を行っています。


最後に緯度も確認しておきましょう。コマンドは書きませんがわかりますよね?


見て明らかなように、北から南にデータが並んでいることがわかります。ちなみに、NCEPや(最近別件でいじった)歴史再解析データは南から北(マジかっ!て思いますよね)。


Monday, 4 November 2024

氷河合宿その5a:NetCDFの中身を確認

 いったんダウンロードから離れて、入手したNetCDFファイルの中身を見ていきましょう。

作業フォルダを作成し、その中に気象データと標高データを入れておきます。練習用&いろいろなコマンドのメモとして、「practice.py」とでも名付けたファイルを作成しましょう。Sublime Textと関連付けしておけば、ダブルクリックですぐに確認、編集できるので便利。

作業フォルダ内でターミナルを立ち上げ、「ipython」と入力&リターンすると、以下のような画面になるはずです。


ここに以下の二行を入力して、NetCDFを使えるようにします。

import numpy as np
import netCDF4

この時点でエラーが出る人は、エラーメッセージを選択コピーして、chatGPTに聞いたりして対処してください。うまくいったら、次のコマンドでダウンロードしたファイルを読み込みます。

nc = netCDF4.Dataset('rw2017mon.nc')

これで、「nc」という配列(行列)の中に、ファイル内のデータが格納されましたので、中身を見ていきます。まずは単純に「nc」とだけ打ってリターン。

nc
バーッと情報が出てきますが、ビビらずに半ばあたりにある「dimensions,,,」以下を選択してコピーし、メモ上にペーストしましょう。パラメータ毎に改行すると、以下のようになっています。

含まれるパラメータ
dimensions(sizes): date(12), latitude(41), longitude(81)
int64 date(date)
float64 latitude(latitude)
float64 longitude(longitude)
float32 u10(date, latitude, longitude):東西風
float32 v10(date, latitude, longitude):南北風
float32 d2m(date, latitude, longitude):露点温度
float32 t2m(date, latitude, longitude):気温
float32 sp(date, latitude, longitude):地表面気圧
float32 tp(date, latitude, longitude):降水量
float32 ssrd(date, latitude, longitude):下向き短波放射
float32 strd(date, latitude, longitude):下向き長波放射

dimensions(sizes): date(12), latitude(41), longitude(81)は、行列の構造です。
u10以下がダウンロードしてきた気象要素で、これがデータを読み出すときに使う名称です。

データ領域
このデータはヒマラヤのロールワリン地域を含む領域で、指定の際に
     30
70 -|- 90
     20
こんな範囲を指定しました。南北10°、東西20°の範囲です。ERA5の解像度は0.25°なので、南北40、東西80のデータがあればよさそうですが、実際には41と81あることがわかります。

これは、データの代表ポイントが間隔0.25°の格子の中心にあるか、四隅に配置されているかの違いで、ERA5は後者(四隅)で、そのために+1こデータがあります。

※最近あまり使ってる論文を見ませんが、NCEPは格子の中央に配置されている。

なので、例えば座標(70,20)のデータが代表する領域は解像度0.25°を半分にした±0.125で、19.875~20.125、69.875~70.125の領域ということになります。これは後々、自分が抽出したい地点がどのグリッドに含まれているか?を求めるのに重要になってきます。

試しにそれぞれのパラメータを確認してみましょう。

nc['u10']

と打ち込むと情報がバーッと出てきますが、ここでチェックしておきたいのは単位(Unit)です。今回入手したパラメータをまとめると

nc['u10'] / nc['v10'] → units: m s**-1
nc['t2m'] / nc['d2m'] → units: K
nc['sp'] → units: Pa
nc['tp'] → units: m
nc['ssrd'] / nc['strd'] → units: J m**-2

Kから℃への換算はよいですね?地表面気圧も普段はhPaを使うので要注意。降水量はmmではなく、mですが、単に1000倍すればいいのか、この段階ではわかりませんね。放射の単位もワットWではなくジュールJですが、一か月分の総量なのでしょうか?確認が必要です。

というわけで、ERA5のOverviewで確認しましょう。

Total Precipitation
選択部分に書かれているように、このデータは一日当たりの降水量をm表記したものです。換算の際は、mmにしたうえで、月ごとの日数をかけてやる必要があります。

Surface Solar Radiation Downwards
こちらも一日当たりの積算エネルギーなので、Wへの換算のためには24*60*60 = 86400で割ってやる必要があります(長波も同じ)。

ずいぶん長くなってしまったので今回はここまで。

Sunday, 3 November 2024

氷河合宿その4b:ERA5標高データ

複数年の自動ダウンロードの前に、「ERA5の標高データ」を入手します。標高の補正に必要なことと、ERA5の座標系がトリッキーなことを確認するために必要です。

どこにあるか?ですが、

前回参照した「Monthly single level」の
「Variable → Other」の中の「Geopotential」

です。「Year」「Month」「Time」はどれか一つだけ選び、APIコードを生成して前回と同じようにダウンロードしておきます。「Geographical area」は「Whole available region」、つまり全球を入手しておきましょう。下記の説明書きにもあるとおり、この要素は時間変化しないので、一度入手しておけば使いまわしができます。ファイル名はお好きなように。

「Elevationじゃないの?」と思った人は良く気づきました。早速、ページの上に戻り、「Overview」のタグを選択し、「Geopotential」で検索してみましょう。

単位とともにどういうパラメータかについての説明がありますね。最後の三行に「地表面標高は重力加速度で割れば得られる。」とあり、重力加速度の値まで提供されています。

つまり、NetCDFに入っているのは標高そのものではないので、9.80665で割る必要がある、ということです。

お次はPythonでNetCDFファイルの中身を見ていくことにします。