※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

PythonでMODISのL3(SMI)を緯度経度で範囲を切り出して図にする。

#モジュールのインストール
from pyhdf.SD import SD,SDC
from matplotlib.pylab import *
from mpl_toolkits.basemap import Basemap

#HDFファイルの読み込み
(フルパスで指定)
#f=SD("C:\\test\A20030602014090.L3m_MC_CHL_chlor_a_4km.hdf", SDC.READ)
f=SD("C:\\test\A20021822002212.L3m_MO_CHL_chlor_a_9km.hdf", SDC.READ)
f.datasets() #アクセス可能なデータを見る
v=f.select("l3m_data") #l3m_dataをvにいれる.

#ここを変更
(#緯度経度)
#換算入力(例:北海道周辺を切り出し)
lah=46   #緯度上
lal=41   #緯度下
lol=139   #経度西
lor=146   #経度東

#ここを変更 #mesh換算計算
(0スタートだからプラスとかいらない)
me=12 #4km=24, 9km=12 #変更
lahv=me*(90-lah)
lalv=me*(90-lal)
lolv=me*(lol+180)
lorv=me*(lor+180)

v3 = v[lahv:lalv,lolv:lorv] #範囲指定

v3.shape #行数と列数の確認

x = linspace(lol,lor,(lor-lol)*me+1) #等間隔のベクトルを作成
y = linspace(lah,lal,(lah-lal)*me+1) #データに比べ一つ多くグリットを作成する

xx,yy = meshgrid(x,y)

#basemapに重ねる
m=Basemap(projection="cyl",
	lat_0=(lah+lal)/2, lon_0=(lol+lor)/2,
	llcrnrlat=lal, llcrnrlon=lol, #lower left corner
	urcrnrlat=lah, urcrnrlon=lor,
	resolution="f", # c l i h f
	area_thresh=1000.0,
	fix_aspect=False)
figure()
m.pcolor(xx, yy, v3, vmin=0, vmax=2) #z軸の調整
m.drawcoastlines()
m.drawparallels(arange(lal,lah,1), labels=[1,0,0,0])
m.drawmeridians(arange(lol,lor,1), labels=[0,0,0,1])
colorbar()
show()
|