Octave 気象データを描画(気温)
やること
アメリカ海洋大気庁(NOAA)/ NCEP の GFS(Global Forecast System)を使って、全球の地表気温を描画します。GFSは、地球全体をカバーした気象予測モデルです。
データ ダウンロード
GFSの予想データをダウンロードします。
今回は、「2021年7月13日 0時(UTC)」を予想したデータをダウンロードしました(ファイル名: gfs_3_20210713_0000_000.grb2)。
コード
grib2形式のファイルをnetcdf形式に変換します。こちらの記事を参考にさせていただきました。ありがとうございました。
NOAAが提供している wgrib2 というツールを使って、grib2をnetcdfへ変換するようです。wgrib2をインストールして、下記のコマンドで変換できます。
# 今回ダウンロードした「gfs_3_20210713_0000_000.grb2」ファイルを、
# netcdf形式に変換して、「gfs_3_20210713_0000_000.nc」というファイル名で出力
wgrib2 gfs_3_20210713_0000_000.grb2 -netcdf gfs_3_20210713_0000_000.nc
netcdfになってしまえば、あとは描画するだけです。
# netcdf を取り扱うパッケージを読み込み
pkg load netcdf
# read data (気温データを読み込み、ケルビンから摂氏へ変換)
in_file = "gfs_3_20210713_0000_000.nc"; #ダウンロードしたGFSデータのPATHを設定
in_info = ncinfo(in_file);
lon = ncread(in_file,"longitude");
lat = ncread(in_file,"latitude");
dat = ncread(in_file,"TMP_surface");
dat=flipud(rot90(dat));
dat = dat - 273.15;
# draw (-90℃~60℃まで1℃間隔で等温線を描画)
contourf(lon,lat,dat,[-90:1:60],'linestyle','none');
colormap(jet);
colorbar('southoutside')
caxis([-5 40]); #カラーバーを -5℃~40℃の範囲にする
daspect([1 1]);
hold on;
# 地形データを読み込み
topo=load("../etopo_1deg.mat"); #地形データのPATHを設定
# 世界地図を描画 (標高0mの等高線)
contour(topo.lon,topo.lat,topo.topo,[0 0],"k");
※世界地図のデータはこちらの記事で作成しています。
描画結果
それでは、お疲れさまでした。