なんとかubuntuにGMTがインストールができた（ホントかな？）ので等高度線を描いてみることにしました
難しいことはわかりませんが、かえってGMTで予想図を描く手順がわかるかもしれません
わかる人が見れば、グチャグチャな手順なんだろうと思います
つまり グチャグチャなのは覚悟してください
GPVデータはターミナル（windowsのDOS窓みたいなやつ）でwgrib2などを使って,データが切り出します
GPVデータのバイナリファイルは
Z__C_RJTD_20150824000000_GSM_GPV_Rgl_FD0000_grib2.bin
を使いました
2015年8月24日００Zイニシャルみたいです
１．バイナリファイルからデータを得る
早崎先生の
によると
wgrib2 ifile.grb2
上記に加え，グリッド数，データの単位，データの範囲(最大・最小値など)を知りたいとき
wgrib2 -V ifile.grb2
さっそく端末（ターミナル）で打ってみる
$ wgrib2 Z__C_RJTD_20150824000000_GSM_GPV_Rgl_FD0000_grib2.bin
結果
1.1:0:d=2015082400:PRMSL:mean
sea level:anl:
1.2:0:d=2015082400:PRES:surface:anl:
1.3:0:d=2015082400:UGRD:10m above ground:anl:
1.4:0:d=2015082400:VGRD:10m above ground:anl:
1.5:0:d=2015082400:TMP:2m above ground:anl:
・・・・
・・・・
1.14:0:d=2015082400:TMP:1000mb:anl:
1.15:0:d=2015082400:VVEL:1000mb:anl:
再びwgrib2使用法によると
２．wgrib2 ifile.grb2 -match ":HGT:" -grib ofile.grb2
ifile.grb2からジオポテンシャル高度だけ切り出して、GRIB-2形式でofile.grb2
に出力．
Z__C_RJTD_20150824000000_GSM_GPV_Rgl_FD0000_grib2.binはうっとうしいので
GPVファイルと書きます
端末で次のように打ち込むと
$ wgrib2 GPVファイル -match":HGT:1000 mb:" -grib 1000HGT.grb2
1000HGT.grb2ができました
さらに1000HGT.grb2の情報をみるため
$ wgrib2 -V 1000HGT.grb2
と入力した結果
1:0:vt=2015082400:1000mb:anl:HGT Geopotential Height [gpm]:
ndata=259920:undef=0:mean=93.4211:min=-403.51:max=441.99
grid_template=0:
lat-lon grid:(720 x 361) units 1e-06 input WE:NS output WE:SN res 48
lat 90.000000 to -90.000000 by 0.500000
lon 0.000000 to 359.500000 by 0.500000 #points=259920
無事できたようですが困りました.
なんとgrb2はGMTでは使えないし変換するソフトもなさそう.
大元のサイト
http://www.cpc.ncep.noaa.gov/products/wesley/wgrib2/netcdf.html
によるとncファイルに変換できそうです.
勿論,私は英語なんてわからいですが運が良かったんですね.
wgrib2 ../example/eta.t00z.awphys18.grb2 -netcdf eta.nc
なんて一文がありました
$ wgrib2 1000HGT.grb2 -netcdf 1000HGT.nc
としてみる 無事にできたが・;
ncはまだGMTでは使えない
grbやgrb2をncに変換してくれるソフトはいくつかあるみたいですが・・ grb2やncをgrbに変換してくれるソフトは見つかりませんでした hoge.ncはGrADSとか言うもので使う形式らしいFORTORUNの世界？
気象屋さんはGMTを使わないのか？
試行錯誤？偶然？で いちからはじめるGMT その１７にあったncdumpを試してみた $ ncdump 1000HGT.nc としてみる.なんと幸いダンプしてくれた.
$ ncdump 1000HGT.nc>1000HGT.txt
1000HGT.txtができた
緯度軽度も入っている0.5刻みのようですが余計（？）なので手作業でカット 面倒ですがtxtファイルでGMTで使えるファイル
（とりあえずやっつけ仕事でmyline2.pyを作成しました）
を作ってGMTでgrdファイルをつくることにする 参考に後で付録してpythonで作ったmyline2.pyを示します 2.GMTで等高度線を描いてみる 予想された等高度線は見ることはできなかったので
いちいち館野や輪島の高層実況を確認していたなあ どの位の高さに０度線があるかも知らないで雪予想なんてできる訳ないのになあと思っていた 1000HGT.txtをとりあえず手作業で修正して1000HGTg.txtをつくり myline2.pyで 1000HGTg.txtをGMTで扱える1000HGTgm.txtへ変換した
次に1000HGTgm.txtを1000HGT.grd（バイナリファイル）に変換します
いちからはじめるGMTその７４をみると
surface signal.txt -Gsignal.grd -R-6/6/-6/6 -I0.125
-Gオプションはgrdファイル名，-Iオプションはグリッド間隔，-Rオプションは緯度経度
元データのグリッドは0.5だからIは半分の0.25にしました領域は日本が入るようにします
$ surface 1000HGTgm.txt -G1000HGT.grd -R120/150/20/50 -I0.25
1000HGT.grdが無事作成されました
それではGMTで書いてみます
図はファイル名1000H.epsにしました $ pscoast -Jm1:5000000 -R135/145/32/38 -W -Bg5 -Gtan -K > 1000H.eps
続けて $ grdcontour 1000HGT.grd -R135/145/32/38 -J -A10 -O >>1000H.eps （ ＞＞1000H.epsであることに注意 ＞1000H.epsではありません）
オプション−Aはコンターの間隔みたいです 無事にできたので gimpでpngに変換（エクスポートして保存）1000H.pngにしました 90度傾いていたので直した図がこれ
バンザイ
北高タイプかな？ついでだから北東気流について少し考えてみよう 850hPaの等高度線がこれ850hPaは関東の沿岸で弱い気圧の谷になってますね やっぱり風が必要か
まだ気象関数は作っていないので味美測候所からお借りします
1000hPaはないので925hPa 925hPa 温度と風たしかに宮城県沖に寒気があって流れこんでいるように見えるな
850hPa 相当温位と風850hPaの前線状のものが関東の925hPa寒期に乗り上げるような形かな？ しっかりみてないのであまりコメントしないほうが・・よいかな？
さらについで 1000hPaの温度線を描いてみる絶対温度になっているようで273度引いて考えてください
本州付近は温度のメリハリがないな 宮城県沖に寒気があって流れこんでくるって説明はどうなんだろう？ 100ｍ位上空の寒気が地上付近に落ちてきたってことになりそうだが・・どうして？ 1000hPaは925hPaに対して不安定？ でも,それは1000hPa温位のほうが高いってことだな？ まあ、カラクリは追々と考えるとして・・
今まで見ることがなかった湿球温度線を見るのが楽しみだな 気象関数はこれから作るので時間がかかります・・ 雪が話題になる頃にはできるかな？ 付録 pythonで作ったmyline2.py
参考にしてください はずかしいけど・・実力がないことがバレバレ
#!/usr/bin/env python #-*- coding: utf-8 -*-
#データを整える関数 def myline(arg):
mylist= arg.split(',')#,でリストに変換 a=[]#新しいリストの初期化 kizami=0.5 latitude=-90 longitude=0 #データ作成のループ for i in mylist: i=i.strip()#空白や改行をなくしている #データの作成 #x座標値 y座標値 １０００hPaの高さ a.append(str(longitude)+' '+str(latitude) +' '+i) #print(str(longitude)+' '+str(latitude)) if longitude >359.5: longitude=0 latitude=latitude+kizami if latitude>89.5: latitude=-90 else: longitude=longitude+kizami #ループは終わり,ここから先は結果をみながら調整してください #リストの最後を表示 print a[-1] #変数のタイプを調べる print type(a) #リストの最後を切る a.pop() #a.strip() #リストから最終的なデータを作る mydata='\n'.join(a) return mydata
#関数終わり
#1000HGTg.txteを読んで関数mylinedeで1000HGTgm.txtを作ります f=open('1000TMP.txt','r') s=f.read() f.close() s=myline(s) f=open('1000TMPgm.txt','w') f.write(s) f.close()
0 件のコメント:
コメントを投稿