地形図を描く

産総研
last updated 2023.9.22
gmt_nav.gif(9191 byte)
作業メモ
GMT3　です
岸本さんの地形ファイルとは：地質調査総合センター研究資料集No.353「海陸を合わせた日本周辺のメッシュ地形データの作成:Japan250m.grd」2000年
起動する
  1. CD-Rをセットする
  2. 扱う有効数字を大きくする
    　　C:\Gmt>gmtset D_FORMAT %f6
  3. ファイルに収録されている範囲を確認する
    　　C:\Gmt>grdinfo d:\japan250m.North.v1a.grd
    　　ここでx_min 128 x_max 148 y_min 30 y_max 46を得る。6400×7680
    　　xinc=0.00312, yinc=0.0020836, z_min -9675, z_max 3720
  4. 同様にSouthもチェックすると
    　　C:\Gmt>grdinfo d:\japan250m.South.v1a.grd
    　　ここでx_min 122 x_max 144 y_min 24 y_max 39.3333を得る。7040×7360
    　　xinc=0.00312, yinc=0.0020836, z_min -9675, z_max 3720
  5. 使う 　
      　
    • ほしいデータを切り出す
      C:\Gmt>grd2xyz d:\japan250m.North.v1a.grd -R144:15/145:15/42:50/43:20 > akkeshi.txt
      範囲の指定の仕方は　西／東／南／北　ですべて度単位の場合は　144.25/145.25/42.83/43.33　となり、度分秒単位のときは上記のようにコロンで区切る。
      こうすると経度（deg）、緯度(deg)、標高（ｍ）のデータができる。
      144.1515636 43.3322926 594.0000006
      144.1546876 43.3322926 546.0000006
      144.1578136 43.3322926 528.0000006
    • 新測地系に変換する
      日本測地系で処理されたファイルを国土地理院で配布しているプログラムTKY2JGD（配布終了）を使って、世界測地系に変換したが、現在は web版が公開されている
前処理
  • 前処理として度単位を度分秒に変換する必要がある。globa_i3.gif(1129 byte)
  • ダウンロードしたプログラムTKY2JGDを起動する。
  • パラメータファイルTKT2JGD.parをセットする
  • 変換形式、緯度経度を選択
  • 一括変換をクリック
  • 別ウインドウが出るので入力ファイル、出力ファイルを指定し、処理開始ボタンをクリック
  • かなり時間がかかるので他の仕事をやって待つ！！
  • できたファイルは以下のような文言がでる（出た）。これは削除する。

    このファイル"akkeshi.gsi.out"は，プログラムTKY2JGD Ver.1.3.79が"akkeshi.gsi.in"を読み込んで計算処理したものです。
    使用した変換パラメータファイルは，"TKY2JGD.par" Ver.2.1.2です。
    次に示すように，各行の最初の２つの数字が，変換されたJGD2000系の緯度，経度を表しています。
    JGD2000系（計算値）
    緯度 経度
    ddmmss.sssss dddmmss.sssss
    行末に「3parameters」があるものは，地域毎のパラメータがなかったか3パラメータで変換するよう設定されていたため，「東京大正」三角点における3パラメータで変換したものです。 また，「-9999.」がある行は，変換されなかった行です。
    以上のどちらでもない行は，「地域毎の変換パラメータ」で変換された行です。
    ただし，コメント行や数値の形式が不正な行は，変換されずにそのまま出力されます。

    残すデータ部分
    432005.27322 1440851.42390
    432005.27345 1440902.67267
    432005.27368 1440913.92144
    430350.51673 1450828.64049
    430350.51711 1450839.88949
    430350.43785 1450850.92162 3parameters
    430350.43826 1450902.17064 3parameters


    　 ここで3parametersの文字も消す。
    エディタを使ってもよいし、awkを使って1項めと2項めだけのファイルを作っても良いー＞akkeshi.gsi2.out
    GMT用に緯度・経度（deg）に変換し、｛経度(deg)、緯度(deg)、標高（ｍ）｝という並びにする。akkeshi.wgs.dataを作成。
    wcコマンドでakkeshi.gsi2.outとakkeshi.txtのデータ数が等しいことを確認する
    データを並び替えて標高値を戻す。globa_i3.gif(1129 byte)

batファイルを作成
  • GMTで図を描くためのbatファイルを作成
    makecpt -Chaxby -T-500/1000/100 -Z > topo.cpt
    blockmean d:\Myhome\yoreyore\akkeshi.wgs.data -R144:15/145:15/42:50/43:20 -I15c -V > akkeshi.blk
    surface akkeshi.blk -R144:15/145:15/42:50/43:20 -I15c -T0.4 -Gakkeshi.grd -V
    grdimage akkeshi.grd -R144:15/145:15/42:50/43:20 -JU55/15 -Ctopo.cpt -P -K -V > akkeshi.eps
    grdcontour akkeshi.grd -R144:15/145:15/42:50/43:20 -JU55/15 -C100 -K -O -V >> akkeshi.eps
    psscale -D5/-1/8/0.2h -Bf50a500g100:topography(m): -Ctopo.cpt -K -O -V >> akkeshi.eps
    pscoast -R144:15/145:15/42:50/43:20 -JU55/15 -Df -W2 -K -O -V >> akkeshi.eps
    psbasemap -R144:15/145:15/42:50/43:20 -JU55/15 -Ba0.5f0.25g0.25WSne -O -V >> akkeshi.eps
メモ

grdcontour：コンター図を描く。読み込むgrdファイル、-Rで描く範囲、 -Jで図法とスケール
-JM15：メルカトール図法で図の横幅1が5cm
-Jm10：メルカトール図法で経度１度あたり10cm
-JU55/15：UTMのゾーン55で図の幅15cm
Zone:52(126E-132E,129E), 53(132E-138E,135E), 54(138E-144E,141E), 55(144E-150E,147E)
-Cでコンター間隔。ここでは100m。以下のようにgrdcontourを2回使ってそれぞれ制限を付けてコンター間隔を変える手もある。

grdcontour akkeshi.grd -R144:15/145:15/42:50/43:20 -JU55/15 -C200 -L200/3000 -P -K -V > akkeshi.eps
grdcontour akkeshi.grd -R144:15/145:15/42:50/43:20 -JU55/15 -C50 -L-1000/0 -K -O -V >> akkeshi.eps

-C200 -L200/3000：コンター間隔200m、200m以上3000ｍ以下について
-C50 -L-1000/200：コンター間隔50m、-1000m以上0ｍ以下について

pscoastの海岸線データとずれてしまうところがあるので標高０ｍはうまく外すのもこつ
-Bで軸
a：annotation：目もり間隔
f：frame：黒白の間隔
g：grid：桝目の間隔

サンプル：コンター間隔10mgloba_i3.gif(1129 byte)