温度の鉛直分布の表示を試みる - その 2 - ライブラリ忘備録編 -

Pygrib Python

　今回は、前回試した気象庁の MSM から一地点のみをわざわざ取り出して鉛直方向に表示するコードの忘備録編です。

　MSM の予報は 5 km メッシュで作成されています。よって、見たい地点の予報データを、その地点の緯度経度を含む 5 km メッシュのデータと見なして、それを取り出すこととしました。

 

Ⅰ. GRIB 2 から必要な範囲のデータを取り出す

　これまで予報データは、水平分布を把握するべく地図上にしか図示しませんでした。そのため、描画範囲は Cartopy の set_extent () で緯度経度を指定していました。つまり、データは全部読み込んでいましたが、描画する範囲のみを限定していた、ということになります。

 

Ⅰ-① Pygrib の範囲指定によるデータ取得

　鉛直方向の描画では Cartopy は使いませんので、 MSM から必要な範囲のみの GRIB2 データを取り出せれば、と思い、Pygrib の説明を読んでみたところ、ありました、ありました。

>>> data, lats, lons = grb.data(lat1=20,lat2=70,lon1=220,lon2=320)

上記の data () を使って、緯度経度を指定、そのデータを取得できます。上記は一例ですが、lat1 = 最小緯度、lat2 = 最大緯度、lon1 = 最小経度、lon2 = 最大経度で囲まれる領域のデータです。

jswhit.github.io

 

Ⅰ-② 気象庁の MSM のメッシュの確認

　では、5 km は緯度経度にすると何度位になるのでしょう。地球 1 周 40000 km として計算でもいいのですが、もっときちんとした数字があるだろうと思い、気象庁の配信資料に関する技術資料を確認します。

www.data.jma.go.jp

　例えば No.500 で確認すると、緯度 0.0625 度、経度 0.05 度とあります。

 

Ⅰ-③ MSM GPV から一メッシュのデータを取り出すプログラム

　ということで、データの格子点が一つは確実に入り、二つは入らないようにするため、緯度経度幅 0.08 度ということで（単にざっくり決めただけですが）、コードを書いてみます。

import pygrib

# Set parameter 
dat = "0923"　　# 9 月 23 日 12Z
tm ="12"
lat=35.764　　　# 見たい場所の緯度経度
lon=140.392
dev=0.04　　　  # 緯度経度幅 0.08 度

lvls=[1000,975,950,925,900,850,800,700,600,500,400,300,250,200,150,100]　　# 気圧面

filenm = "Z__C_RJTD_2021"+dat+tm+"0000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin"
grbs = pygrib.open(filenm)

def get_value(req_data):
    val = grbs.select(forecastTime=0, level=lvl)[req_data]
    val_pos = val.data(lat1=lat-dev, lat2=lat+dev, lon1=lon-dev, lon2=lon+dev)
    Val_pos=float(val_pos[0])
    val_set.append(Val_pos)

val_set=[]
for lvl in lvls:    
    req_data = 3　　# 温度
    get_value(req_data)

print(val_set)

　上記は MSM GPV の気圧面データを利用するものなので、地上面や他の GPV ではファイル名やデータの種類はちょっと変わると思います。さて、上記のプログラムを動かしてみると以下のようになりました。無事 16 の基準気圧面の温度予報値が取得できています。

[295.95611572265625, 296.2651062011719, 296.1998291015625, 297.4837951660156, 295.5482177734375, 291.1206359863281, 286.7886962890625, 280.7910461425781, 274.5665283203125, 266.0440979003906, 255.32080078125, 239.46371459960938, 228.833740234375, 217.5061492919922, 211.73158264160156, 205.81198120117188]

　ひとまず目標は達成されました。そして、これは個人的希望ですが、鉛直方向の温度分布を見るのであれば、気温の他に露点が欲しいと思いました。もっとも、MSM GPV に気温データはあるのですが、露点データはありません。ですが、相対湿度データはあるので、ここから、露点を求めたいと思います。

 

Ⅱ. 露点を計算する（と言っても、ライブラリを使うだけです）

　普段目にする、例えば空港気象情報には、気温と露点が並んで表示されているので、気温を測るように、露点も何らかの手段で難なく測れるものと思っていたのですが、検索をかけてみると、そう簡単というわけではないようです。事実、何らかの計算や、理論式からの換算表などが用いられることが多いようです。そこで、今回理屈の理解はひとまず断念し（笑）、その値を求める方法のみを確認することとします。　そして、こういう時こそ、python ではライブラリの出番です。

 

Ⅱ-① Metpy により露点を求める

　前回書いたように、Metpy というライブラリは気象関連の計算描画がいろいろ出来るライブラリなので、その簡易版をわざわざ自分で作るのは、本来であれば、単にプログラミングの練習程度にしかなりません。ですが、今のところは描画の段階でなぜかエラーが出るので、今回は、その計算機能部分だけ拝借することとします。

　こういうライブラリであれば当然なのでしょうが、気温と相対湿度から露点を計算する機能があるのはさすが素晴らしい限りです。

unidata.github.io

　Metpy の面白い所は、物理量に対する姿勢が厳格というか、物理量に次元（単位）を付けて計算するところです。実にもっともな話ではあるのですが、電卓の延長のように計算しようとすると、単位分だけコードが長く、ややくどい感は拭えません。おそらく、次元をコードに組み込めた方が利点のある人がユーザーの多くを占めるのでしょうね。

　次元をコードに組み込むことについて書かれたチュートリアルもバッチリあります。

unidata.github.io

 

Ⅱ-② 計算してみましょう

　では、あるまだ穏やかめな暑い夏の日の露点を求めてみましょう。

　気温約 32 ℃、相対湿度 85 ％ の時の露点を求めるコードは以下の様になります。

from metpy.units import units
import metpy.calc as mpcalc

dew_pos = mpcalc.dewpoint_from_relative_humidity(305 * units.kelvin, 0.85)
dewp=dew_pos/units.kelvin -273.15

print(dewp, type(dewp))

　上の結果は以下の通りです。

29.01570497604706 dimensionless <class 'pint.quantity.build_quantity_class.<locals>.Quantity'>

　露点約 29 ℃。湿度の暑い夏は、数度気温が下がったところで、大気がすぐ飽和するので、ベタベタ感がなくなることがないことが分かります。日中多少気温が上がったところで、秋の空気はやはり素敵ですね。

　コード中、units.kelvin を掛けたり割ったりしている点が、先に触れた通り、温度という物理量の次元部分です。露点 29.015.... は無次元（dimensionless）、数値だけですが、長いクラス名が付けられています。細かいことは分かりませんが（笑）、Ⅰ-③ で求めた温度のリスト val_set と同等に扱えます。

 

　前回の記事を書いた時には、コードの忘備録を残そうと意気込んでいたのですが、少し時間が経っただけで、何を残しておこうと思ったのか、そもそも何で忘備録を作ろうと思ったのかすら見事に怪しくなりました。

　あとは取り出した GRIB2 データを図示するだけなので、描画ライブラリに詳しい方の方が色々とできるかと思いますが、自分でも使い勝手の良い図をちょっと模索中なので、また後日、続きのまとめを書くかもしれません。

 

　まとめをほとんど作らない人生を送っている身には、忘備録を作ることは、やはりなかなか骨が折れます。

　本日も最後までお付き合いどうもありがとうございました。