QGISとGDALでラスタ画像からタイル地図を作成する方法ひととおりメモ #GIS

はじめに

スキャンした古地図を良い感じにタイル地図のデータに変換する必要があったので、その都度調べたメモを残します。GISソフトウェアを扱うのはほぼ初めてです。同じ境遇の方の参考になればと思います。

なお、QGISを使ってGCPsを設定する部分以外は、なるべく下記のコマンドを使って処理するようにしています。

gdal_translate
gdalwarp
gdal2tiles.py

システム環境

システム環境は下記のとおりです。

Mac OS X 10.11.3
QGIS 2.14.3-Essen: GISソフトウェア

QGISは公式サイトからダウンロードのリンクを進んでいくとMac OS X installers for QGISのリンクに辿り着きます。するとインストーラ QGIS-*.dmg があるのでこれでインストールできます。

このパッケージに必要なプログラム一式が入っているので、GDAL、NumPy、matplotlib、QGISという説明通りの順番でインストールを進めていきます。QGISの機能の多くは、裏方でPythonやGDALなどを実行していますので、QGISだけでなく、他のツール群も必要になるわけです。

地理参照の概要

下記サイトによると、スキャンした画像を地理情報として扱えるようにするための作業を ジオリファレンシング （地理参照）と呼ぶそうです。

スキャナーを使って取り込んだ地図は、例えばjpegやtiffなどの形式で画像として保存されますが、そのままでは地理情報ではなく画像情報です。画像情報を地理情報にするための作業が地理参照（ジオリファレンシング）であり、QGISはそのための機能を提供しています。
GCNオープンソースGIS勉強会 › QGISによる自然環境情報解析 › 第1章 › デジタイジング：ラスタデータの取り込みと幾何補正

初心者がQGISを使ったジオリファレンシングを行うにあたっては、上記サイトを一読すると理解が早いと思います。私の場合は、ジオリファレンシングという用語を知るのに結構な時間がかかってしまい、最終的にはこのサイトを初めに読むのが良かったと思います。

地図の表示

画像を地理情報として扱えるようにするには、画像のピクセル座標と地図の座標の対応付けをする必要があります。このような、画像の座標と地図の座標の対応を表現する点群を GCP (Ground Control Point; グラウンドコントロールポイント) と呼びます。

今回は、元の画像には地理情報が方角以外は書き込まれていなかったので、QGISを使って画像と地図を見比べながらGCPsを地道に打っていくことにします。そのためには、まず画像との比較対象にする地図を表示する必要があります。

地図を表示するのに下記の方法を試してみました。

TileLayer Plugin: WMTS (Web Map Tile Service) 形式で公開されている地図を表示できます。地図タイルと呼ばれているようです。タイル画像がキャッシュされるので使っているうちに動作がキビキビします。
OpenLayers Plugin: OpenLayersというWMTS形式の地図を表示するJavaScriptのライブラリを使ったプラグインです。上記と同様に表示できますがレスポンスが微妙でした。ただし、Google MapsなどはWMTS形式のタイル画像のURLに直接アクセスすることは禁止されているので上記のTilelayer Pluginは使用できませんが、こっちは違う方法を使ってGoogle Mapsの地図も表示できます。

今回は 地理院地図 と OpenStreetMap を使うことにします。理由は下記の通りです。

今回は城や神社といった代表的な点だけの対応付けで必要十分
となると城や神社が掲載されていて利用が簡単でレスポンスの良い地図が良い
地理院地図のレスポンスはけっこう良い
Google MapsをOpenlayers Pluginで表示するとレスポンスが微妙
OpenStreetMapのサーバはレスポンスが微妙でしたが今回の地図の範囲ではかなり詳細に情報が書き込まれていたのでこちらも利用
TileLayer Pluginはタイル地図のデータのキャッシュが効くので動作レスポンスが良好

TileLayer Pluginをインストールするには、QJISのメニューバーから、プラグイン、プラグインの管理とインストール、TileLayer Pluginを検索、プラグインをインストールを選べばOKです。

TileLayer Pluginで扱う地図は設定ファイルに設定を記述することで追加できます。作者があらかじめ地理院地図の設定ファイルを下記URLで公開しておられますので、利用規約を読んでからTileLayer Pluginの設定ファイルが置かれるディレクトリにコピーすればOKです。

GitHub Gist minorua/GSIMaps.tsv

Mac OS Xの場合は下記コマンドでダウンロードしてコピーできます。

$ wget https://gist.githubusercontent.com/minorua/7654132/raw/16ab51b5b469995f28a2fcf10b02e35ac04aa9fa/GSIMaps.tsv -P ~/.qgis2/python/plugins/TileLayerPlugin/layers/

デフォルトでは、XYZFrameとTMSFrameというタイル地図の枠だけ表示するレイヤーしかありませんが、上記のファイルを追加することで、地理院地図が使えるようになります。

さらに下記の行を追加しておけばOpenStreetMapのタイル地図も使えるようになります。

OpenStreetMap   OpenStreetMap   https://a.tile.openstreetmap.org/{z}/{x}/{y}.png        1       1       18

地図を表示するには、QJISのメニューバーから、Web、タイルレイヤプラグイン、タイルレイヤを追加する…、から、例えば、GISMaps、標準地図などを選択して追加ボタンを押せば良いでしょう。

以上で地図が無事表示できましたはずです。

ジオリファレンシング

GDALジオリファレンサープラグインの起動

QGISには「GDALジオリファレンサー」というプラグインが標準で入っていました。まずはこれを有効化します。

GDALジオリファレンサー Pluginを有効化するには、QGISのメニューバーから、プラグイン、インストール済み、GDALジオリファレンサーのチェックボックスをONにして閉じます。Mac OS Xの場合、チェックボックスの色が薄くて本当にONになっているか不安ですが、色が薄いことに特別な意味は無くちゃんとONになっているようです。（たまに、QGISを起動し直した場合に、プラグインが無効化されていることがあります。これはなかなか面倒です。）

GDALジオリファレンサープラグインが有効化されていれば、QGISのメニューバーから、ラスタ、ジオリファレンサー、ジオリファレンサーを選択します。すると新しいウィンドウが開きます。このウィンドウに画像を読み込んで、元のウィンドウに表示されている地図と画像との対応関係 (GCPs) を設定していくわけです。

画像の読み込み

ツールバーの一番左にラスタを開くボタンがありますのでこれでラスタ画像を指定します。

（とりあえず練習用に自分が住んでいる近くのラスタ画像を探してみようと思ったのですがパブリックドメインの画像はなかなかありませんでした。観光地図とかオープンデータで転がっているかなと思いましたが見つかりませんでした。ジオパークの地図なんて良かったのですが無断転載はダメと記載があったのでやめました。）という理由があって、とりあえず今回は、テキサス大学が公開している1945-1946年に米軍が作成した日本の都市計画地図を使ってみます。なんと鳥取県の鳥取市も掲載されていました。（ちなみに、この地図画像を正確にTMS (Tile Map Service) 形式にデータ化しようとしている方がおられました（テキサス大学図書館米軍のJapan City Plans TMSデータ公開）。鳥取はまだみたいなので、ここまでの精度を出す元気はないので、自分で適用にやってみます。

空間参照系の選択

ラスタ画像を読み込もうとすると「空間参照システム選択」「レイヤRaterのCRSを指定してください」というダイアログがでてきます。つらいです。なんのことやという気持ちです。ということで調べてみました。QGISによる自然環境情報解析の空間参照系の解説が解りやすかったので必要な部分だけを以下に整理します。

空間参照系: 空間参照システムと呼ぶこともあります。英語では SRS (Spatial Reference System) と呼びます。
コード体系: 空間参照系はいくつかの構成要素で定義されるためそれらの組み合わせを識別するためのコードが定義されています。
- EPSGコード: EPSG (European Petroleum Survey Group) という団体が整理している。
- ESRIコード: ESRI (Environmental Systems Research Institute) というArcGISの開発元が整理している。
地理座標系: 英語では CRS (Coordinate Reference System) と呼びます。

ここで選択するEPSGコードは、QGISでの地図表示に使っているEPSGコードに揃えて選びます。QGISの右下のステータスバーに表示されています。ここで違うEPSGコードを選ぶと、ジオリファレンサーの設定にある「完了時にQGISにロード」した時に表示がずれてしまいます。と言っても読み込んだあとで変更できるので大した問題ではありません。ここでは EPSG:3857 を選びました。

画像と地図のリンク

まず画像を読込んだら変換の設定のアイコンをクリックして出力ラスタのファイル名を指定します。ファイル名を指定しておくと、画像と地図のGCPsを設定するときに、画像と地図の場所、ズームレベルを双方向にリンクさせることができるようになり、GCPsの設定が楽になります。リンクはジオリファレンサープラグインのツールボタンから設定できます。傾きはリンクされないみたいなので、地図の傾き（QGISの左下あたり）を微調整すると楽です。

幾何補正

変換の設定で選べる変換タイプによって必要なGCPsの数が違います。必要な点数を決めると誤差が表示されます。取り急ぎ誤差を見たい場合は必要な点数が少なく3点で済む多項式1変換を選ぶと良いと思います。

ちなみに私がQGISで試してみたところ下記のようになりました。なお、処理の都合上、コマンドラインを使いたいので、GDALでは未サポートの変換タイプについては調べていません。カッコ内は gdalwarp に与えられるコマンドラインオプションです。

1 線形変換: GDALでは未サポート
2 ヘルマート: GDALでは未サポート
3 多項式1変換: ( -order 1 -co): 3点
4 多項式2変換: (-order 2 -co): 6点
5 多項式3変換: (-order 2 -co): 10点
6 シンプレートスプライン: (-tps -co):
7 投影変換: 4点 GDALでは未サポート

これらの変換タイプの概要についてはArcMapのラスターデータセットのジオリファレンスの基礎にまとまっていました。

また、下の番号ほど良い複雑な変換になるそうです。複雑な変換だからといって良い結果になるとは限らないとのことです。（デジタイジング：ラスタデータの取り込みと幾何補正）

QGISは表記揺れが気になりますが、それぞれの意味は下記のとおりでしょうか。

変換タイプ: 幾何補正
多項式1変換: 一次アファイン多項式
多項式2変換: 二次アファイン多項式
多項式3変換: 三次アファイン多項式

今回は、画像を単純にシフト、サイズ変更、回転するだけで十分なので 多項式1変換（一次アファイン多項式） を使うことにしました。

リサンプリング方法

下記の種類がありました。リサンプリング方法と表記されていますが、アルゴリズムは内挿法（ないそうほう、Interpolation）や補間法とも言うようです。

最近傍
線形
キュービック
キュービックスプライン
ランチョシュ

例のごとく、略語が使われているので正確な名称と用途を下記URLの1、2、3から調べたのをまとめます。

QGIS表記	日本語	gdalwarpの `-r` オプション	用途
最近傍	最近隣内挿法	near (default)	最も高速なリサンプリング手法。各値がクラス、メンバー、または分類 (土地利用、土壌タイプ、または森林タイプなどのカテゴリデータ) を表す名目データまたは順序データに使用してください。ごつごつした出力。（1より）
線形	共一次内挿法	biinear	最近隣内挿法を使用する場合よりも滑らかなサーフェスが得られます。連続的なサーフェスに適しています。標高、傾斜角、空港からの騒音公害、河口付近における地下水の塩分濃度などはすべて、連続サーフェスとして表現できる事象であり、共一次内挿法を使ってリサンプリングするのが最適です。滑らかな結果。（1より）
キュービック	三次たたみ込み内挿法	cubic	共一次内挿法よりもデータが鮮明になる傾向にあります。したがって、このリサンプリング手法は、航空写真や衛星画像などの画像のリサンプリングによく使用されます。最も鮮明な結果。（1より）
キュービックスプライン	3次スプライン	cubicspline	画像効果をスムーズにしますがぼかし効果が強くなる場合があります。（2より）
ランチョシュ	ランチョス法（ランツォシュ法）	lanczos	最もシャープな画像に仕上がりますが、不自然な結果に仕上がることもあります。（2より）

今回は実は古地図のデータなので キュービック を使おうと思います。

GCPsの追加

あとはGDALジオリファレンサーと地図を行ったり来たりしながらGCPを打っていくだけです。横に並べると作業がしやすいです。まだ、誤差が表示されるので、あまりに誤差が大きい点は、採用しない方が良いと思います。イラスト地図などのように測量的な正確性を求めない場合は別です。

それと、画像の4つの角をGCPとして追加しておくと良いと思います。あとで、画像を回転させたときに、角が切れてしまう場合があります。普通は gdalwarp コマンドのオプション -ts で指定できるとは思うのですが、うまくいかなかったので、これで対処しています。

ジオリファレンスの開始

変換の設定で、変換先SRSを EPSG3857にして、完了時にQGISにロードするにチェックをいれて、ジオリファレンスの開始をクリックします。うまくいけば、地図にオーバーレイされるはずです。

背景が黒い場合は、レイヤパネルからラスタデータのレイヤを右クリックをして透過性の透過ピクセルリストで色をピックアップすれば透明になります。

画像が出ない場合は、SRSの指定を間違えている可能性があるので、その場合は、変換先SRSを変えてもう一度同じ処理をしてみるか、上記と同じように右クリックして、一般情報の空間参照システムの指定を変えてみると良いと思います。

GDALスクリプトの生成

GDALジオリファレンサーのGDALスクリプトの生成というボタンを押すと下記のスクリプトが生成されます。

$ gdal_translate -of GTiff -gcp 4642.52 2440.49 1.49423e+07 4.20057e+06 -gcp 6077.98 2385.98 1.49423e+07 4.20067e+06 -gcp 5996.21 568.949 1.49422e+07 4.20076e+06 -gcp 4676.59 423.587 1.49421e+07 4.20068e+06 -gcp 58.4274 13.4271 1.49419e+07 4.20037e+06 -gcp 14.4113 7026.66 1.49424e+07 4.20001e+06 -gcp 9918.04 7026.66 1.49429e+07 4.2007e+06 -gcp 9852.01 13.4271 1.49424e+07 4.20106e+06 "/Users/masayuki/Documents/git.higashino.center.tottori-u.ac.jp/chizu-chizu/chizu-chizu/images/map-orig.jpg" "/var/folders/cc/470s1lj92dx_7z9ykv36kt0m0000gn/T/map-orig.jpg"
$ gdalwarp -r cubic -order 1 -co COMPRESS=NONE  "/var/folders/cc/470s1lj92dx_7z9ykv36kt0m0000gn/T/map-orig.jpg" "/Users/masayuki/Documents/git.higashino.center.tottori-u.ac.jp/chizu-chizu/chizu-chizu/images/map-orig_modifiedx.tif"

gdal_translate: ラスタデータ異なるラスタデータに変換するコマンドです。 -gcp を指定することで、GCP情報付きのラスタデータに変換します。
gdalwarp: 画像の投影法と移動を行うユーティリティコマンドです。GCP情報付きのデータをもとに、幾何補正とリサンプリングを行えます。

GDALスクリプトの修正

前述したスクリプトを少し修正します。

`gdal_translate`

まず、-of オプションでGTiffと書いているのに出力ファイル名はjpgなので違和感があります。拡張子はjpgですが実際は -of オプションの指定どおり、ファイル形式はGeoTIFF形式です。また、ファイルサイズが 3.8MBから 199MBと約52倍になっており、GeoTIFF形式である必要性がない場合は、ディスク容量がかなり無駄です。

そこでVRT形式を使います。VRT形式の拡張子は .vrtで、中身はXMLファイルです。コマンドは下記のとおりです。

$ time gdal_translate -of VRT -gcp 4642.52 2440.49 1.49423e+07 4.20057e+06 -gcp 6077.98 2385.98 1.49423e+07 4.20067e+06 -gcp 5996.21 568.949 1.49422e+07 4.20076e+06 -gcp 4676.59 423.587 1.49421e+07 4.20068e+06 -gcp 58.4274 13.4271 1.49419e+07 4.20037e+06 -gcp 14.4113 7026.66 1.49424e+07 4.20001e+06 -gcp 9918.04 7026.66 1.49429e+07 4.2007e+06 -gcp 9852.01 13.4271 1.49424e+07 4.20106e+06 ./map1.jpg map1.vrt
Input file size is 9921, 7016

real    0m0.031s
user    0m0.006s
sys 0m0.020s

結構長いですが下記のようなXMLファイルが生成されます。オリジナルの画像ファイルのファイル名（ファイルパスでもOK）が含まれているので、ファイルを置く場所には注意が必要です。

map.vrt

<VRTDataset rasterXSize="9921" rasterYSize="7016">
  <Metadata>
    <MDI key="EXIF_ColorSpace">1</MDI>
    <MDI key="EXIF_Orientation">1</MDI>
    <MDI key="EXIF_PixelXDimension">9921</MDI>
    <MDI key="EXIF_PixelYDimension">7016</MDI>
    <MDI key="EXIF_ResolutionUnit">2</MDI>
    <MDI key="EXIF_XResolution">(600)</MDI>
    <MDI key="EXIF_YResolution">(600)</MDI>
  </Metadata>
  <Metadata domain="IMAGE_STRUCTURE">
    <MDI key="INTERLEAVE">PIXEL</MDI>
  </Metadata>
  <GCPList>
    <GCP Id="" Pixel="4642.5200" Line="2440.4900" X="1.494230000000E+07" Y="4.200570000000E+06" />
    <GCP Id="" Pixel="6077.9800" Line="2385.9800" X="1.494230000000E+07" Y="4.200670000000E+06" />
    <GCP Id="" Pixel="5996.2100" Line="568.9490" X="1.494220000000E+07" Y="4.200760000000E+06" />
    <GCP Id="" Pixel="4676.5900" Line="423.5870" X="1.494210000000E+07" Y="4.200680000000E+06" />
    <GCP Id="" Pixel="58.4274" Line="13.4271" X="1.494190000000E+07" Y="4.200370000000E+06" />
    <GCP Id="" Pixel="14.4113" Line="7026.6600" X="1.494240000000E+07" Y="4.200010000000E+06" />
    <GCP Id="" Pixel="9918.0400" Line="7026.6600" X="1.494290000000E+07" Y="4.200700000000E+06" />
    <GCP Id="" Pixel="9852.0100" Line="13.4271" X="1.494240000000E+07" Y="4.201060000000E+06" />
  </GCPList>
  <VRTRasterBand dataType="Byte" band="1">
    <Metadata />
    <ColorInterp>Red</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">./map1.jpg</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="9921" RasterYSize="7016" DataType="Byte" BlockXSize="9921" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="9921" ySize="7016" />
      <DstRect xOff="0" yOff="0" xSize="9921" ySize="7016" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="2">
    <Metadata />
    <ColorInterp>Green</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">./map1.jpg</SourceFilename>
      <SourceBand>2</SourceBand>
      <SourceProperties RasterXSize="9921" RasterYSize="7016" DataType="Byte" BlockXSize="9921" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="9921" ySize="7016" />
      <DstRect xOff="0" yOff="0" xSize="9921" ySize="7016" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="3">
    <Metadata />
    <ColorInterp>Blue</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">./map1.jpg</SourceFilename>
      <SourceBand>3</SourceBand>
      <SourceProperties RasterXSize="9921" RasterYSize="7016" DataType="Byte" BlockXSize="9921" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="9921" ySize="7016" />
      <DstRect xOff="0" yOff="0" xSize="9921" ySize="7016" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

`gdalwarp`

画像が回転されてデータがない部分を透過処理したいので -dstalpha オプションを加えます。出力ファイル形式はGeoTIFFです。拡張子は .tif です。

$ gdalwarp -r cubic -order 1 -dstalpha -co "COMPRESS=NONE" map1.vrt map1.tif
Creating output file that is 11556P x 11950L.
Processing input file map1.vrt.
0...10...20...30...40...50...60...70...80...90...100 - done.

各リサンプリング方法と処理時間と出力サイズは下記のとおりでした。出力ファイルサイズは同じなんですね。画質は、うーん、あんまり変わらない気がしますが、拡大した時に near よりも cublic の方が滑らかな気がします。

Resampling methods	Time	Output File Size
near	0m27.782s	527MB
bilinear	0m36.178s	527MB
cubic	0m46.902s	527MB
cubicspline	1m28.051s	527MB
lanzos	1m15.282s	527MB
average	0m36.785s	527MB
mode	0m41.318s	527MB

`gdal2tiles.py`

下記コマンドを実行すると gdal モジュールが参照できないようです。

$ gdal2tiles.py -h
Traceback (most recent call last):
  File "/Library/Frameworks/GDAL.framework/Programs/gdal2tiles.py", line 46, in <module>
    import gdal
ImportError: No module named gdal

次のコマンドによるとpython2.7を使用しているようです。

$ head -n 2 `which gdal2tiles.py`
#!/usr/bin/env python2.7
# -*- coding: utf-8 -*-

私の環境ではMacPortsでインストールしたpython2.7 (/opt/local/bin/python2.7) も入ってそっちのPATHを優先していたので、QGISのインストーラを信じれば、おそらく次の /usr/bin/python2.7 を使えばよさそうです。

$ which -a python2.7
/opt/local/bin/python2.7
/usr/bin/python2.7

ということでこうしました。普通はwhichではなくフルパスを指定した方がよいと思います。

/usr/bin/python2.7 `which gdal2tiles.py` -h

結局、こういうコマンドになります。ソースファイルのSRSは EPSG:3857 なので --s_srs epsg:3857 オプションを指定しています。

time /usr/bin/python2.7 `which gdal2tiles.py` --s_srs epsg:3857 map1.tif
Generating Base Tiles:
0...10...20...30...40...50...60...70...80...90...100 - done.
Generating Overview Tiles:
0...10...20...30...40...50...60...70...80...90...100 - done.

real    0m36.816s
user    0m34.957s
sys 0m1.189s

ちなみに地図タイルの総ファイルサイズは 43MB でした。

おわりに

まとめると下記のスクリプトで画像からタイル地図が生成できました。GCPsの値さえあれば、ほぼ自動化できそうです。

$ gdal_translate -of VRT -gcp 4642.52 2440.49 1.49423e+07 4.20057e+06 -gcp 6077.98 2385.98 1.49423e+07 4.20067e+06 -gcp 5996.21 568.949 1.49422e+07 4.20076e+06 -gcp 4676.59 423.587 1.49421e+07 4.20068e+06 -gcp 58.4274 13.4271 1.49419e+07 4.20037e+06 -gcp 14.4113 7026.66 1.49424e+07 4.20001e+06 -gcp 9918.04 7026.66 1.49429e+07 4.2007e+06 -gcp 9852.01 13.4271 1.49424e+07 4.20106e+06 ./map1.jpg map1.vrt
$ gdalwarp -r cubic -order 1 -dstalpha -co "COMPRESS=NONE" map1.vrt map1.tif
$ /usr/bin/python2.7 `which gdal2tiles.py` --s_srs epsg:3857 map1.tif