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

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

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



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

  

久しぶりに綺麗なグラフとかを描きたくて
GMTを触ってみた

バージョンが５に上がっていて、いろいろ
変わったところが多い。
インストールはubuntuの場合cmakeでしたり
するようになっていたり。

プログラミングもコマンドからモジュールという
概念に変わっているようだ。
例えば
gmt pscoast ほにゃらら
と入力する。

より日本語の情報が少ないので
できたことはアップしてみる。

今回は紀伊半島周辺の海岸線図に
県境と、距離スケールと方位を示す矢印を追記する。
矢印、というかベクトルを書くのに手間取った。

gmtのバージョン5以降は
インストールできているとして、
県境のデータをダウンロード。
こちらを参考に
http://www.geocities.jp/ne_o_t/otenki.htm

後は、観測地点なんかも記載できたらいいので、
std_name.dataと名付けたファイル
　　NachiRiver 33.668 135.904
と、std＿point.dataと名付けたファイル
　　135.904 33.668 LM 0.0 12p,Helvetica,255/0/0 NachiRiver
を用意。

今回はmacで作業してるので、あんまりgmtのファイルが
どこにあるかわからないので、shareと名付けたファイルに
県境のデータ ken.txtを置いて、以下のようなシェルスクリプトを書いて実行。

#!/bin/bash
# kenzakai_wakayama.sh

rangeSW=134.5/136.5/33/35
scale=M10c
afg=a1f0.5g0
ruler=136/33.3/33.2/50k
pen1=0.6p,black #[width,color]
pen2=0.4p,50
g_color=200
s_color=255

point_data=stn_point.dat
name_data=stn_name.dat

fig=Wakayama.eps

gmt gmtset FORMAT_GEO_MAP +DF
gmt set MAP_VECTOR_SHAPE 0

gmt psbasemap -J$scale -R$rangeSW -B$afg -X3c -Y3c -P -V -K > $fig
gmt pscoast -J -R -B$afg -Lf145/30/35/500 -Dh -P -W$pen1 -G$g_color -S$s_color -L$ruler -V -O -K >> $fig
gmt psxy ../share/ken.txt -J -R -W$pen2 -O -K >> $fig
gmt psxy $point_data -i2,1 -J -R -Sc0.2c -G255/0/0 -V -O -K >> $fig #地点の位置
gmt pstext $name_data -J -R$rangeSW -D0.5/0 -F+j+a+f -O -K >> $fig #地点の名前 -Dはシフト

#方向矢印の記述
#矢印
gmt psxy -R -J -SV0.4+s+e+a45 -W1 -G0 -O -K << EOF >> $fig #片矢印
135,33.3,135,33.65 #始点、終点
EOF
#横棒
gmt psxy -R -J -SV0.5+s -W0.8 -G0 -O -K << EOF >> $fig #線分を引く
134.93,33.45,135.07,33.45 #始点、終点
EOF
#Nの文字
gmt pstext -R$rangeSW -J -D0/0.1 -F+j+a+f -O << EOF >> $fig
135 33.65 BC 0.0 12p,Helvetica,black N
EOF

#白抜きのpngの作成
#convert -density 100 test2.ps test2.png
# mogrify -trim -density 300 -bordercolor 'transparent' -border 10x10 -format png *.eps
#mogrify -trim -density 300x300 -transparent white -format png $fig


以上。
最後でコメントアウトしてるけど、
imagemagicを使って、白のところを透明にしたpngファイルも
作成できるようにしている。

こんな図ができると思う。（クリックして拡大）


ベクトルを書くのに、psxyの-SVオプションが結構戸惑った。
WANtaroHP (Programing on Mac)さん
http://kbtkt.web.fc2.com/sub_0_gmt5.html
を参考にさせてもらった。ありがとうございます。  

礫の移動を表示するため
QGISに、
WKTのテキストデータを
レイヤー追加して、ラインを引いた。

矢印表示にできるように
ラインのプロパティを以下のように設定。



こんな感じです。（赤の矢印）


画像をクリックしたら、大きく表示されます。  

３年ほど使ってるMac mini(Late2012)の
ハードディスクをSSDに交換してみました。
ずっと思ってたけど、ちょっと難しく思っていたので。

このために買ったのは、

と、
交換で出てくるHDDのための外付けケース

と、
ちょっと特殊なドライバーが必要だということで


SSDが値上がり気味ってことなので、
慌てて買いました。

件のMacmini

一年ほど前にメモリは８GBに交換済み。

裏返すと

黒い蓋をヒネって開ける。


購入したドライバーでビスを外して、
ごそごそとファンを外す。




金属の金網みたいものを外す。この時
wifiのアンテナ線も外します。


すると、黒いカバーに覆われたHDDが見える。
ごそごそと引き抜く。
自分のMacminiは、下側（ちゃんと置くと）にあったので、取りやすい方だった。
上側にあるものもあるらしい。そうなると、作業がさらに増えるみたい。



HDDについてる部品を外して、SSDに装着して、
手順を遡って、元に戻すだけ。


でも、wifiのアンテナの接続部分とか、
fanの電源端子が小さくて、つけるのは結構手間取った。
ピンセットがあると便利。

SSDには予めSierraを外付けのHDDケースを使って、
インストールおいたので、すぐに起動できた。
MacはOSを簡単にダウンロードして、外付け
のドライブにインストールできるのがすごいな。

この交換作業は色々ネットに情報があるので、
今更ですが。
Mac miniはLate2012に根強い人気があるらしい。