はじめに

今回は備忘録として日本の地形の3D画像を作成する手順について書く．

環境

OS windows 10
cygwin 2.11.0-1
(cygwinのパッケージとして，Perl，shells，GDALをインストールしておく)
GMT 4

基盤地図情報ビューア（ページ下部）

やり方（備忘録）

まず作りたい範囲の地形のメッシュデータを国土地理院の基盤地図情報ダウンロードサービスからダウンロード．
（この際にアカウント登録等の手続きが必要）
ダウンロード後，基盤地図情報ビューアを開き，ダウンロードしたメッシュデータをzipファイル形式のまま，ドラッグアンドドロップ．
すると以下のように地形データが可視化される．

その後，画面上部，エクスポート/標高メッシュをシェープファイルへ出力/　からシェープファイルを出力．この際当地域が直角座標系の何系にあたるか調べておく．
出力された4ファイルの位置で，cygwinから以下のコマンドを実行する．

ogr2ogr -f GMT ***_gcs.gmt ***.shp

上記コマンドからバイナリ形式のシェープファイルをGMTで利用可能なテキスト形式に変換．***にはファイル名．

出力されたファイルをgmtというフォルダに格納，同じ階層で以下のソースコードを実行する．

gmt2grd.pl

#!/usr/bin/perl

# 球面またはデカルト座標系のxyz形式データをgrdファイルに変換する

$ratio = 2.0; # 勾配データ出力倍率

$input_path = $ARGV[0];

@paths = split (/\//, $input_path);

$gmt_dir = join ("/", @paths[0 .. $#paths - 1]);

$filename = $paths[$#paths];

if ($filename =~ /(.+)_gcs(.*)\.gmt$/)
{
    # 地理座標系の場合
    $filename = $1 . "_gcs" . $2;
#   $dx = 2.25 / 3600.0; # 50m
    $dx = 2.25 / 10 / 3600.0; # 5m
#   $dy = 1.5 / 3600.0; # 50m
    $dy = 1.5 / 10 / 3600.0; # 5m
    $prec = "7f";   # 出力データにおける小数点以下の精度
    $r = 0.1 * $dx;
}
elsif ($filename =~ /(.+)_xy(.*)\.gmt$/)
{
    # 平面直角座標系の場合
    $filename = $1 . "_xy" . $2;
    $dx = "5";  # 出力データの格子幅（単位：メートル；x方向）
    $dy = "5";  # 出力データの格子幅（単位：メートル；y方向）
    $prec = "2f";   # 出力データにおける小数点以下の精度
    $r = $ratio;
}
else
{
    print "Usage: gmt2grd.pl [ data_name ]_[ gcs | xy* ]\r\n";
    exit;
}

上記ソースコードを

perl gmt2grd.pl gmt/***_gcs.gmt

で実行．
実行後，gmt/grdの中に***_gcs.grdファイルが作成されているのを確認．
次に同じ階層でgrdフォルダを作成し，作成したgrdファイルを格納．そして以下のソースコードをgmtフォルダと同階層で実行．

grd2eps.sh

#!/usr/bin/csh

# 3D作図用のtcshスクリプト

touch ./grd/null.xyz

xyz2grd ./grd/null.xyz -G./grd/null.grd -R./grd/###_gcs.grd -V -N-2

makecpt -Crelief -T-301/299/10 -Z -D > dem.cpt

grdmath ./grd/***_gcs.grd ./grd/null.grd AND = tmp10.grd
grdgradient tmp10.grd -D -Stmp10_slope.grd
grdmath tmp10_slope.grd 0.00001 MUL NEG = tmp10_slope.grd
#以下の-Rオプションは該当地域の地理座標系緯度経度を西/東/南/北の順に入力
#gmtフォルダに出力された***_gcs.txtに緯度経度が書き込まれている
grdview tmp10.grd -Jm190c -Jz0.0025c -
R*********/********/********/******** -Cdem.cpt -E90/30 \
    -Itmp10_slope.grd -Qc600 -N0/127/127/127 -Wf0.25p -Bg2000/g2000WeSn -V -Xc -Yc -P \
    > ***_gcs_3D_1.eps

上記ソースコード内のgrdviewコマンド -Rオプションに地理座標系の緯度経度を入力する必要があるので注意．
grdviewコマンドの各オプションを設定することで，3D画像の表示方法を自分で変えられる．-Jzでは鉛直方向の縮尺，-Eでは方位角，仰角を設定．他詳細はこちら．

このようなeps形式の3D画像ができれば完成．