和歌山サイエンスカフェ「インフィニティ」

GMT5.3.1を使い、国土地理院の基盤地図情報サービスの
１０mメッシュのDEMから、地形の陰影図や鳥瞰図を作る。

まずはデータのダウンロード
http://www.gsi.go.jp/kiban/
（今回は紀伊半島の一部のデータを利用）

登録してダウンロードし解凍してできる
FG-GML-5035-37-dem10b-20161001.xml
を利用します。

GMTで扱えるgrd形式に変換する必要がある。
grd形式に変換するスクリプトを
産業技術研究所の野田さんが公開されているので
利用させていただく。ありがとうございます。
https://staff.aist.go.jp/a.noda/programs/jpgis2grd/jpgis2grd.html

このスクリプトを動かすとgrd形式に変換できるが
一部環境に合わせて変更する。
（commandを「gmt xyz2grd」とか出来上がるファイル名などの変数を）

$ perl jpgis2grd.pl
を実行
しかし、エラーで止まる。
途中のどこかに３列目の数値がないところがあるようだ。
（今回は 843751行目）

このスクリプトのglobal optionに
my $outfile2 = "out2.xyz"
を追記。

grdに変換するあたりを

print STDOUT "Making a grd file.\n";
# my $command="gmt xyz2grd $outfile -R$xmin/$xmax/$ymin/$ymax -I$xint/$yint -G$grdfile -N-9999 -Dlongitude/latitude/height/1/0/$opts{title}/$opts{remark} ";
my $command_sed="sed -e '843751d' $outfile > $outfile2";
my $command="gmt xyz2grd $outfile2 -R$xmin/$xmax/$ymin/$ymax -I$xint/$yint -G$grdfile -di-9999 -Dlongitude/latitude/height/1/0/$opts{title}/$opts{remark} ";
if ($opts{verbose}) { print STDERR "Running \"$command\"\n"; };
system("$command_sed");
system("$command");
system("gmt grdinfo $grdfile");

と一部修正。
（perlに詳しくないので、あんまりいいやり方ではなさそうですが）

再度、
$ perl jpgis2grd.pl -v -n
を実行
（実行状況や中間ファイルを残すのに　-v -nオプションをつけた）

すると、out.grdファイルができるので
とりあえず、これを
FG-GML-5035-37-dem10b-20161001.grd
と名前を変更。

次はgmtで陰影図を描く
以下のスクリプトで
#!/bin/bash
#
# DEM10mメッシュのgrdファイルから陰影図を作る。
# 2016/12/12
#
# perl jpgis2grd.pl で出力される範囲をコピー
# xmin=135.87500000, xmax=136.00000000, ymin=33.58333333, ymax=33.66666667

grd_file=FG-GML-5035-37-dem10b-20161001.grd
shade_file=5035-37.shade
eps_file=5035-37_inei.eps
cpt_file=color.cpt
ctb=-1000/1000/50 #color table min/max/interval
region=135.875/136.00/33.5833333/33.6666667

gmt makecpt -Cglobe -T$ctb -Z > $cpt_file

gmt gmtset FORMAT_GEO_MAP D

gmt grdgradient $grd_file -A90 -G$shade_file -Ne0.6 -V
gmt grdimage $grd_file -R$region -JM15 -C$cpt_file -I$shade_file -X4 -Y4 -K > $eps_file
gmt psscale -D7.5/-1/12/0.3h -C$cpt_file -Ba200f100:"altitude(m)": -I -V -K -O >> $eps_file
gmt psbasemap -R$region -JM15 -Ba0.05f0.025g0.01 -O -V >> $eps_file

こんな図ができる



次は鳥瞰図
#!/bin/bash
#
# DEM10mメッシュのgrdファイルから陰影図を作る。
# 2016/12/12
#
# perl jpgis2grd.pl で出力される範囲をコピー
# xmin=135.87500000, xmax=136.00000000, ymin=33.58333333, ymax=33.66666667

grd_file=FG-GML-5035-37-dem10b-20161001.grd
shade_file=5035-37.shade
d3_file=5035-37_3d.eps
cpt_file=color.cpt
ctb=-1000/1000/50 #color table min/max/interval
region=135.875/136.00/33.5833333/33.6666667
boundary_cap=wSnEZ
gmt makecpt -Cglobe -T$ctb -Z > $cpt_file
gmt gmtset FORMAT_GEO_MAP D
gmt grdgradient $grd_file -A90 -G$shade_file -Ne0.6 -V
gmt grdview $grd_file -Jm130 -Jz0.002 -C$cpt_file -I$shade_file -p150/25 -Ba0.05f0.025g0.01 -Bza400 -B$boundary_cap -Qi300 -N0+ggray -V > $d3_file

を実行すると、鳥瞰図ができる。



角度とかは適宜設定する。