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