2015年8月26日水曜日

GMTで等高度線を描いてみました


なんとかubuntuGMTがインストールができた(ホントかな?)ので等高度線を描いてみることにしました
難しいことはわかりませんが、かえってGMTで予想図を描く手順がわかるかもしれません
わかる人が見れば、グチャグチャな手順なんだろうと思います
つまり グチャグチャなのは覚悟してください


GPVデータはターミナル(windowsDOS窓みたいなやつ)でwgrib2などを使って,データが切り出します

GPVデータのバイナリファイルは
Z__C_RJTD_20150824000000_GSM_GPV_Rgl_FD0000_grib2.bin
を使いました
2015824日00Zイニシャルみたいです


1.バイナリファイルからデータを得る

早崎先生の

によると


1.変数名,日付などの基本情報だけを知りたいとき
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使用法によると
2.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
無事できたようですが困りました.
なんとgrb2GMTでは使えないし変換するソフトもなさそう.
大元のサイト
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では使えない
grbgrb2ncに変換してくれるソフトはいくつかあるみたいですが・・
grb2ncgrbに変換してくれるソフトは見つかりませんでした
hoge.ncGrADSとか言うもので使う形式らしいFORTORUNの世界?
 
気象屋さんはGMTを使わないのか?
 
試行錯誤?偶然?で
いちからはじめるGMT その17にあったncdumpを試してみた
$ ncdump 1000HGT.nc
としてみる.なんと幸いダンプしてくれた.
 
$ ncdump 1000HGT.nc>1000HGT.txt
 
1000HGT.txtができた
緯度軽度も入っている0.5刻みのようですが余計(?)なので手作業でカット
面倒ですがtxtファイルでGMTで使えるファイル
(とりあえずやっつけ仕事でmyline2.pyを作成しました)
を作ってGMTgrdファイルをつくることにする
参考に後で付録してpythonで作ったmyline2.pyを示します




2.GMTで等高度線を描いてみる
予想された等高度線は見ることはできなかったので
いちいち館野や輪島の高層実況を確認していたなあ
どの位の高さに0度線があるかも知らないで雪予想なんてできる訳ないのになあと思っていた

1000HGT.txtをとりあえず手作業で修正して1000HGTg.txtをつくり
myline2.pyで 1000HGTg.txtGMTで扱える1000HGTgm.txtへ変換した
 
次に1000HGTgm.txt1000HGT.grd(バイナリファイル)に変換します
 
いちからはじめるGMTその74をみると
 
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はコンターの間隔みたいです
無事にできたので
gimppngに変換(エクスポートして保存)1000H.pngにしました
90度傾いていたので直した図がこれ
 
バンザイ
 
北高タイプかなついでだから北東気流について少し考えてみよう
850hPaの等高度線がこれ




















850hPaは関東の沿岸で弱い気圧の谷になってますね
やっぱり風が必要か 
まだ気象関数は作っていないので味美測候所からお借りします
 
1000hPaはないので925hPa
925hPa 温度と風























たしかに宮城県沖に寒気があって流れこんでいるように見えるな
 850hPa 相当温位と風























850hPaの前線状のものが関東の925hPa寒期に乗り上げるような形かな?
しっかりみてないのであまりコメントしないほうが・・よいかな? 
さらについで
1000hPaの温度線を描いてみる絶対温度になっているようで273度引いて考えてください
 
 
 
 
 
 
 
 
 
 
 
 
 
本州付近は温度のメリハリがないな
宮城県沖に寒気があって流れこんでくるって説明はどうなんだろう?
100m位上空の寒気が地上付近に落ちてきたってことになりそうだが・・どうして?
1000hPa925hPaに対して不安定?
でも,それは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座標値 1000hPaの高さ 
 
      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を読んで関数mylinede1000HGTgm.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 件のコメント: