地図タイルをGMTの地図に貼り込む

加納靖之

2021年3月14日 16:33

前回，気象庁の地震月報の震源データから，2016年熊本地震の震央の時空間分布図をつくってみました．

GMTの海岸線データを使ったわけですが，地理院地図などのタイル地図に重ねて描くこともできます．タイル地図の画像をダウンロードして貼り込むだけです．

もちろんGISソフトやWeb GISを使って描画できますが，GMTはGMTの良さがあるし，GMTのほうが慣れているという方もあると思います．なお縮尺を大きくするときは，座標系や震源決定の精度などいろいろご注意ください．

対象とするデータ

気象庁「地震月報(カタログ編)」 震源データ
地理院タイル 淡色地図
OpenStreetMapタイル OpenStreetMap Foundation Japan 「osm.jpタイルサーバの切り替えについて」

使ったツール・環境

gmt6 https://www.generic-mapping-tools.org/
gdal_translate (GDAL 3.0.4, released 2020/01/28)
awk
WSL2上のubuntu Ubuntu 20.04.2 LTS 「Windows 10 用 Windows Subsystem for Linux のインストールガイド」

描き方

震源データの取得

まず気象庁のホームページの震源データのページから震源データをダウンロードしましょう．今回は2016年熊本地震の震央分布を書きたいので2016年のデータ（h2016）です．解凍しておきます．

地図タイルの取得

地図タイルをダウンロードします．あとで描き体地図の範囲をもとに，必要な地図タイルを取得します．そのためのスクリプトは下のようになります．ここで，座標（緯度・経度）と地図タイル番号の換算はOpenStreetMap Wikiの「Slippy map tilenames」という記事を参考にしました．

この記事にシェルスクリプト＋awkでの処理方法も書いてあります．環境によってawkの小数点以下の処理が違う（私の環境では切り捨て）部分を修正しています．また，日本周辺だけ扱うときはあまり気にならないと思いますが，マイナスの地図番号にならないように修正しました．

#!/bin/sh
# 地図タイルの操作のための関数
# https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
# を修正（負の緯度経度，awkの小数点以下の扱い）
xtile2long()
{
  xtile=$1
  zoom=$2
  echo "${xtile} ${zoom}" | awk '{printf("%.9f", $1 / 2.0^$2 * 360.0 - 180)}'
}

long2xtile() 
{ 
  long=$1
  zoom=$2
  echo "${long} ${zoom}" | awk '{ $1+=$1<0?360:0;
  xtile = ($1 + 180.0) / 360 * 2.0^$2; 
  #xtile+=xtile<0?-0.5:0.5;
  xtile+=xtile<0?-1.0:0.0;
  printf("%d", xtile ) }'
}

ytile2lat()
{
  ytile=$1;
  zoom=$2;
  tms=$3;
  if [ ! -z "${tms}" ]
  then
    # from tms_numbering into osm_numbering
    ytile=`echo "${ytile}" ${zoom} | awk '{printf("%d\n",((2.0^$2)-1)-$1)}'`;
  fi
  lat=`echo "${ytile} ${zoom}" | awk -v PI=3.14159265358979323846 '{ 
  num_tiles = PI - 2.0 * PI * $1 / 2.0^$2;
  printf("%.9f", 180.0 / PI * atan2(0.5 * (exp(num_tiles) - exp(-num_tiles)),1)); }'`;
  echo "${lat}";
}

lat2ytile() 
{ 
  lat=$1;
  zoom=$2;
  tms=$3;
  ytile=`echo "${lat} ${zoom}" | awk -v PI=3.14159265358979323846 '{ 
  tan_x=sin($1 * PI / 180.0)/cos($1 * PI / 180.0);
  ytile = (1 - log(tan_x + 1/cos($1 * PI/ 180))/PI)/2 * 2.0^$2; 
  #ytile+=ytile<0?-0.5:0.5;
  ytile+=ytile<0?-1.0:0.0;
  printf("%d", ytile ) }'`;
  if [ ! -z "${tms}" ]
  then
  # from oms_numbering into tms_numbering
    ytile=`echo "${ytile}" ${zoom} | awk '{printf("%d\n",((2.0^$2)-1)-$1)}'`;
  fi
  echo "${ytile}";
}

# タイルの取得
## 引数
ZOOM=$1
LONGMIN=$2
LONGMAX=$3
LATMIN=$4
LATMAX=$5

## 取得するタイル番号の上限と下限
XTILEMIN=`long2xtile $LONGMIN $ZOOM`
XTILEMAX=`long2xtile $LONGMAX $ZOOM`
YTILEMIN=`lat2ytile $LATMIN $ZOOM`
YTILEMAX=`lat2ytile $LATMAX $ZOOM`
echo $XTILEMIN $XTILEMAX $YTILEMIN $YTILEMAX

# 地理院タイル（淡色地図）
BASEURL='https://cyberjapandata.gsi.go.jp/xyz/pale'
# OpenStreetMap Japan
#BASEURL='https://tile.openstreetmap.jp'
for I in `seq $XTILEMIN $XTILEMAX`;
do
 for J in `seq $YTILEMAX $YTILEMIN`;
 do
   ## 地理院地図からダウンロード
   URL=$BASEURL/$ZOOM/$I/$J.png
   echo $URL 
   curl --silent $URL -o ${ZOOM}_${I}_${J}.png
   ULX=`xtile2long $I $ZOOM`
   ULY=`ytile2lat $J $ZOOM`
   I1=`echo $I | awk '{print $1+1}'`
   LRX=`xtile2long $I1 $ZOOM`
   J1=`echo $J | awk '{print $1+1}'`
   LRY=`ytile2lat $J1 $ZOOM`
   ## 座標をつけてGeoTIFFに変換
   gdal_translate -a_ullr $ULX $ULY $LRX $LRY ${ZOOM}_${I}_${J}.png ${ZOOM}_${I}_${J}.tif
 done
done

引数として，ズームレベル，左端経度，右端経度，下端緯度，上端緯度を与えます．ズームレベルはとりあえず9でどうでしょう．ズームレベル7で九州の4/5くらいが入るので，それより2段階ズームしたもの，という感じです．ズームレベルと描画範囲の関係は，国土地理院の「タイル座標確認ページ」が参考になると思います．

/bin/sh gettiles.sh 9 130 132 32 33.5

地理院タイルはpng画像です．gdal_translateを使ってダウンロードしたタイル画像に座標の情報を入れ，GeoTIFF形式に変換して保存しています．ズームレベル_Xタイル番号_Yタイル番号.tifというファイルができていると思います．画像ビューアでも確認できます．

地図を描く

さてデータが準備できたので，GMTで地図を書きましょう．スクリプトは以下のようになります．前回のスクリプトに地図タイルの座標とタイル番号の変換用の関数と，タイルを貼り込む部分を追加しています．まさにタイルを貼るように，順番に貼りつけています．

#!/bin/sh
# 地図タイルの操作のための関数
# https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
# を修正（負の緯度経度，awkの小数点以下の扱い）
xtile2long()
{
  xtile=$1
  zoom=$2
  echo "${xtile} ${zoom}" | awk '{printf("%.9f", $1 / 2.0^$2 * 360.0 - 180)}'
}

long2xtile() 
{ 
  long=$1
  zoom=$2
  echo "${long} ${zoom}" | awk '{ $1+=$1<0?360:0;
  xtile = ($1 + 180.0) / 360 * 2.0^$2; 
  #xtile+=xtile<0?-0.5:0.5;
  xtile+=xtile<0?-1.0:0.0;
  printf("%d", xtile ) }'
}

ytile2lat()
{
  ytile=$1;
  zoom=$2;
  tms=$3;
  if [ ! -z "${tms}" ]
  then
    # from tms_numbering into osm_numbering
    ytile=`echo "${ytile}" ${zoom} | awk '{printf("%d\n",((2.0^$2)-1)-$1)}'`;
  fi
  lat=`echo "${ytile} ${zoom}" | awk -v PI=3.14159265358979323846 '{ 
  num_tiles = PI - 2.0 * PI * $1 / 2.0^$2;
  printf("%.9f", 180.0 / PI * atan2(0.5 * (exp(num_tiles) - exp(-num_tiles)),1)); }'`;
  echo "${lat}";
}

lat2ytile() 
{ 
  lat=$1;
  zoom=$2;
  tms=$3;
  ytile=`echo "${lat} ${zoom}" | awk -v PI=3.14159265358979323846 '{ 
  tan_x=sin($1 * PI / 180.0)/cos($1 * PI / 180.0);
  ytile = (1 - log(tan_x + 1/cos($1 * PI/ 180))/PI)/2 * 2.0^$2; 
  #ytile+=ytile<0?-0.5:0.5;
  ytile+=ytile<0?-1.0:0.0;
  printf("%d", ytile ) }'`;
  if [ ! -z "${tms}" ]
  then
    # from oms_numbering into tms_numbering
    ytile=`echo "${ytile}" ${zoom} | awk '{printf("%d\n",((2.0^$2)-1)-$1)}'`;
  fi
  echo "${ytile}";
}

# 地震月報の震源ファイルをWIN final形式に変換する関数
jma2final () {
  awk '{
    printf "%04d %02d %02d %02d %02d %5.2f %7.4f %8.4f %6.2f %3.1f %s\n",\
    substr($0,2,4),substr($0,6,2),substr($0,8,2),\
    substr($0,10,2),substr($0,12,2),substr($0,14,2)+substr($0,16,2)/100,\
    substr($0,22,3)+(substr($0,25,2)+substr($0,27,2)/100)/60,\
    substr($0,33,4)+(substr($0,37,2)+substr($0,39,2)/100)/60,
    substr($0,45,3)+substr($0,48,2)/100,
    substr($0,53,2)/10,
    substr($0,69,24);
  }'
}

# 図のサイズ
MAPWIDTH=12
TIMEWIDTH=14

# 描画範囲
YMIN=32.0
YMAX=33.5
XMIN=130
XMAX=132
TICK=a1f0.25
DRY=white
WET=grey80

# WIN final形式に変換，必要な範囲のデータだけ抽出，一時ファイルに保存
cat h2016 | \
jma2final | \
awk '$2==4{print $0}' | \
awk '$8>='$XMIN'&&$8<='$XMAX'&&$7>='$YMIN'&&$7<='$YMAX'{print $0}' > tmp.$$
gmt begin kumamoto_tile
gmt set MAP_FRAME_TYPE plain MAP_DEGREE_SYMBOL degree \
FORMAT_DATE_MAP mm/dd PS_MEDIA a4 \
FONT_ANNOT_PRIMARY 12p,Helvetica,black \
FONT_LABEL 16p,GothicBBB-Medium-EUC-H,black

# 震央分布
gmt basemap -JM$MAPWIDTH -R$XMIN/$XMAX/$YMIN/$YMAX -B$TICK -X3 -Y15

## 地図タイルを貼る
ZOOM=9
XTILEMIN=`long2xtile $XMIN $ZOOM`
XTILEMAX=`long2xtile $XMAX $ZOOM`
YTILEMIN=`lat2ytile $YMIN $ZOOM`
YTILEMAX=`lat2ytile $YMAX $ZOOM`

for I in `seq $XTILEMIN $XTILEMAX`;
do
  for J in `seq $YTILEMAX $YTILEMIN`;
  do
    gmt grdimage ${ZOOM}_${I}_${J}.tif
  done
done

## 線分A-B（時系列を投影する）
gmt plot <<END -W
131.4 33.4
130.4 32.2
END

## M<5
cat tmp.$$ | \
awk '$10<5{print $0}' | \
awk '$7>='$YMIN'&&$7<='$YMAX'&&$8>='$XMIN'&&$8<='$XMAX'&&$9<=20{print $8,$7,$10/(10-$10)/5}' | \
gmt plot -Sc -Wthinnest,grey50

## M>=5
cat tmp.$$ | \
awk '$10>=5{print $0}' | \
awk '$7>='$YMIN'&&$7<='$YMAX'&&$8>='$XMIN'&&$8<='$XMAX'&&$9<=20{print $8,$7,$10/(10-$10)/5}' | \
gmt plot -Sc -Wthinner,red

## ラベルA-B
gmt text << END -F+f12p,0,black+jLT
131.4 33.4 A
130.4 32.2 B
END

## 凡例
gmt plot <<END -Sc -Wthinner,black -N
132.2 32.7 0.467
132.2 32.5 0.3
132.2 32.3 0.2
END

gmt text << END -F+f12p,0,black+jCM -N
132.2 32.62 M7.0
132.2 32.42 M6.0
132.2 32.22 M5.0
END

# 地理院タイル
gmt text << END -F+f10p,GothicBBB-Medium-EUC-H,black+jLT -N
130.75 31.87 地理院タイルに地震月報による震央分布を追記
END

# OSMタイル
#gmt text << END -F+f10p,0,black+jCT -N
#131 31.87 Base map and data from OpenStreetMap and OpenStreetMap Foundation
#END

# 時系列
## 線分A-Bの長さ
DIST=`gmt mapproject -G131.4/33.4+uk <<END | awk '{print $3}'
130.4 32.2
END
`
gmt basemap -JX$TIMEWIDTH/-8 -R2016-04-14T/2016-05-01T/0/$DIST -Bxafg+l"2016年" -Byafg+l"距離, km" -BWSne \
-X-1 -Y-10

## M<5
cat tmp.$$ | \
awk '$10<5{print $0}' | \
awk '$7>='$YMIN'&&$7<='$YMAX'&&$8>='$XMIN'&&$8<='$XMAX'&&$9<=20 {print $8,$7,$10/(10-$10)/5,$1"-"$2"-"$3"T"$4":"$5":"$6}' | \
gmt project -C131.4/33.4 -E130.4/32.2 -W-25/25 -Lw -Q | \
awk '{print $4,$5,$3}' | \
gmt plot -Sc -Wthinnest,grey50

## M>=5
cat tmp.$$ | \
awk '$10>=5{print $0}' | \
awk '$7>='$YMIN'&&$7<='$YMAX'&&$8>='$XMIN'&&$8<='$XMAX'&&$9<=20 {print $8,$7,$10/(10-$10)/5,$1"-"$2"-"$3"T"$4":"$5":"$6}' | \
gmt project -C131.4/33.4 -E130.4/32.2 -W-25/25 -Lw -Q | \
awk '{print $4,$5,$3}' | \
gmt plot -Sc -Wthinner,red

## ラベルAとB
gmt text << END -F+f12p,0,black+jLT -N
2016-5-1T12 0 A
2016-5-1T12 $DIST B
END
gmt end

# 一時ファイルを削除
rm tmp.$$

できあがりはこんな感じです．自分で描けば，もっと高解像度の画像が得られます．下の地図は地理院タイル（淡色地図）に気象庁地震月報（カタログ編）の震源データを重ねたものです．

OpenStreetMapのタイルを貼り込む

上記の地図タイルをダウンロードするスクリプトでは，OpenStreetMap Foundation Japanから地図タイルをダウンロードすることもできます．BASEURLという変数のコメントアウトの部分を変えてください．ここを修正すれば，「標準地図」など他の地理院タイルや，同様の形式で提供されている他のタイル地図も使えます．タイル画像の保存名は同じなので，そのまま実行すると，地理院地図のタイル画像のファイルは上書きされます．