GNU Octave で世界地図を描画
はじめに
気象・海洋・気候変動などの分野で国際的に広く使われているデータフォーマットであるnetCDFを、GNU Octaveというソフトウェアを使って、簡単に可視化させることを目的としています。なおGNU Octaveはオープンソースであり、GPLの条件下で無料で使用することができます。
今回は、地形データをダウンロードして、世界地図を描画させていきます。プログラミングの学習に時間をかけず、簡単にデータを可視化させることに重点を置いています。機械学習や高度な解析が必要であれば、次のステップとして、Pythonなどを学習するとよいでしょう。
インストール
ここでは windows10を使います。
GNU Octave のインストールについては、こちらの記事を参照してください。
標高データの準備
ETOPO1地形データを使用して世界地図を描画させます。
ETOPO1はNOAA（アメリカ海洋大気庁）で公開されている陸地・海底の標高データです。緯度経度1分間隔の2次元格子データとなっています。
下記のページからダウンロードしてみましょう。
ETOPO1ダウンロードページ
今回は南極とグリーンランドの氷床下の地球表面のデータを選択
ETOPO1_Bed_g_gmt4.grd.gz をダウンロード
gzip圧縮してあるので解凍を行い、解凍後の「ETOPO1_Bed_g_gmt4.grd」ファイルを、作業フォルダに移動しておきましょう。ここでは作業フォルダを「 D:\octave\ 」とします（※作業フォルダは好みの場所に作成してください）。
世界地図を描画
Octaveを起動して「現在のディレクトリ」を、作業フォルダにして、ファイルブラウザに「ETOPO1_Bed_g_gmt4.grd」が表示されることを確認しましょう。
以下のコマンドをOctaveのコマンドウィンドウに入力してエンターを押しましょう。
# netcdf を取り扱うパッケージを読み込み
pkg load netcdf
# ETOPO1読み込み
topo_file = "ETOPO1_Bed_g_gmt4.grd";
x = ncread(topo_file,"x");
y = ncread(topo_file,"y");
z = ncread(topo_file,"z");
z=flipud(rot90(z));
[xm,ym] = meshgrid(x,y);
# 1分格子(21601x10801)を1度格子(360x180)に変換
lon=[-179.5:1:179.5];
lat=[-89.5:1:89.5];
[xi,yi]=meshgrid(lon,lat);
zi = interp2(xm,ym,z,xi,yi);
#西経180°から始まるデータを 東経0°から始まるデータに順番を入れ替える
zi=[zi(:,lon>=0) zi(:,lon<0)];
lon=[lon(lon>=0) lon(lon<0)];
lon(lon<0)=lon(lon<0)+360;
topo=zi;
# 描画 (標高0mの等高線)
contour(lon,lat,topo,[0 0]);
daspect([1 1]);【補足】#記号からはじまる行は、コメントと呼ばれるメモで、プログラムとして認識されない行です。
Figure1ウィンドウが新たに開き、世界地図を描画したと思います。
描画された世界地図
【補足】1度格子の地形データを再利用するのであれば、このデータを保存するのが便利です。
# 作成した1度格子の地形データを etopo_1deg.mat というファイル名で保存
save("etopo_1deg.mat","lon","lat","topo");