2016年12月12日

GMT5.3.1で陰影図や鳥瞰図を作る

GMT5.3.1を使い、国土地理院の基盤地図情報サービスの
10mメッシュの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
を実行
しかし、エラーで止まる。
途中のどこかに3列目の数値がないところがあるようだ。
(今回は 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

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



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



タグ :gmt5.3.1grdmac

同じカテゴリー(スタッフ日記)の記事画像
GMT5.3.1で矢印とか地図をかく
QGISで矢印を表現させる
Mac mini Late2012のHDDをSSDに換装
zenfone3 の購入
河川シミュレーションソフトiRIC
土壌雨量指数grib2形式の読み込みプログラム
同じカテゴリー(スタッフ日記)の記事
 2017初頭のウルトラブック比較 (2017-02-26 19:08)
 Latexをwindowsにインストール (2017-02-18 21:28)
 Javaを使った有限体積法による数値計算(1) (2017-02-16 21:07)
 GMT5.3.1で矢印とか地図をかく (2016-12-09 00:08)
 QGISで矢印を表現させる (2016-12-07 20:38)
 Mac mini Late2012のHDDをSSDに換装 (2016-12-06 23:26)

Posted by 和歌山サイエンスカフェインフィニティ at 15:55│Comments(0)スタッフ日記
上の画像に書かれている文字を入力して下さい
 
<ご注意>
書き込まれた内容は公開され、ブログの持ち主だけが削除できます。