おしらせ

オープンソースGISでベクタデータの幾何補正

2014年07月11日

オープンソースGISでベクタデータの幾何補正を行う手順。

  • 主な機会は、形状のゆがみ、位置のズレがある時
  • 主な対象は、ライン、ポリゴン

ちなみに、ラスタデータの幾何補正を行う方法は、検索で多数ヒットする。

  • QGISの「ジオリファレンサー」プラグインを用いた方法
  • GDALコマンド(gdal_translate -gcp、gdalwar -s_srs)を用いた方法
  • GRASS 6.4.3の「幾何補正処理」を用いた方法

利用するツール

OSGeo4Wとそれに含まれる、GDAL/OGRライブラリ、QGISなど
http://trac.osgeo.org/osgeo4w/wiki/OSGeo4W_jp
OSGeo4Wは、アップデートも容易で、多数のパッケージが用意されている。
オープンソースGISを始めるのに便利なツール。

GDAL/OGRライブラリで「GDAL 1.10.0」から利用できるようになった機能がある。
http://www.gdal.org/ogr2ogr.html
それが、ogr2ogr -gcp
現状「GDAL 1.11.0」になっているが、機能は変わっていない模様。

予め取得したGCPの値を用いて、以下のようなコマンド文で幾何補正が可能である。

ogr2ogr ^
-gcp もと座標X もと座標Y 補正座標X 補正座標Y ^
-gcp もと座標X もと座標Y 補正座標X 補正座標Y ^
-gcp もと座標X もと座標Y 補正座標X 補正座標Y ^
-gcp もと座標X もと座標Y 補正座標X 補正座標Y ^
-f "ESRI Shapefile" hosei.shp moto.shp

次のオプションで変換タイプが指定できる模様。
-order n　1次多項式～3次多項式
-tps　シンプレートスプライン

利用するデータ

室蘭市のむろらんオープンデータライブラリ
http://www.city.muroran.lg.jp/main/org2260/odlib.php

「都市計画現況図平成23年版」、「津波浸水深さ」の２つを利用した。
「都市計画現況図平成23年版」は基図情報、「津波浸水深さ」は主題情報と仮定し選定した。
サンプルデータとして軽くするために、室蘭駅周辺にクリップした。
補正対象として、もとの世界測地系から、日本測地系に変換したデータを用意した。
以下の画像は、用意したデータをQGIS2.2.0で開いたところ。
黒色の基図、青色の主題が、もとの世界測地系のデータ。
「＊_tky」のピンク色の基図、オレンジ色の主題が、日本測地系に変換したデータ。

用意したデータ

QGISのオンザフライCRS変換の設定

QGISで異なる空間参照システム（CRS）のデータを開くと、「オンザフライCRS変換」が有効になっているため、最初に開いた方に合わせて表示される。
上記の画像は、「オンザフライCRS変換」を無効にした結果。
「オンザフライCRS変換」を無効にする手順は、

    1. QGIS右下にある地球アイコンをクリック
      CRSステータス
    2. プロジェクトのプロパティ｜CRS　画面が開く
      上部にある、「オンザフライCRS変換を有効にする」のチェックをはずす
      プロジェクトのプロパティ｜CRS

ちなみに、空間参照システム（CRS）の一覧から任意のCRSを選択すれば、QGIS上で投影するCRSとして設定できる。

補足

EPSGコード：世界中の空間参照システムを定義しているコード。
EPSGコードを現在管理しているイギリスのOGPのサイトで検索できる。
http://www.epsg-registry.org/
ほかにも
http://www.spatialreference.org/
QGISの上記の画像にあるプロジェクトのプロパティ｜CRS　画面でも、「フィルター」欄にキーワードを入れれば、該当するEPSGコードが表示される。

幾何補正の手順その１

GCPを、QGISの「座標キャプチャ」プラグインで取得する。
    1. ベクタデータを追加する
      追加した結果は、利用するデータで紹介した画像のようになる
      ベクタデータを追加する

    2. 「座標キャプチャ」プラグインを有効にする
      プラグインを有効にする
      座標キャプチャプラグイン
    3. 必要に応じてCRSを設定する
      座標キャプチャ画面が現れ、２つの空欄がある
      上欄は、任意に表示する座標を変更できる。デフォルトではWGS84（EPSGコード：4326）が設定されている
      下欄は、プロジェクトのプロパティ｜CRS　画面で設定されている座標値が表示される
      CRSの設定
    4. 「キャプチャ開始」
      「キャプチャ開始」ボタンをクリックすると、マウスカーソル上の座標値を随時取得、赤枠に表示されるようになる
      表示された座標値は、「クリップボードへコピー」ボタンでコピー、テキストエディタなどに記録しておく
      キャプチャ開始
GCP取得のポイント

ラスタデータの幾何補正と同様。
対象となるベクタデータの全体をバランスよく囲えるように取得する。
・上下左右の対角線上の交差点中心部などで４点分
・加えて、斜めの対角線上の交差点中心部などで４点分
など、ベクタデータの状態により柔軟に対応する。
４点でも十分な場合もある。多くのGCPを取ることに注力しない方が良い。
基準点のポイントデータ（座標値）がある時は同時に読み込み、
基準点の位置（座標値）をGCPとして取得した方が、絶対座標値の向上につながる。

GCPを取得した結果

以下は、補正対象となるピンク色の基図上に赤丸の４点GCPを取得したところ。

補正するGCPの全体箇所補正するGCPの左上補正するGCPの右上
補正するGCPの右下補正するGCPの左下

以下は、補正先となる黒色の基図上に赤丸の４点GCPを取得したところ。もとが同じ基図なので、同じ場所の取得が容易であった。

補正先のGCPの全体箇所

補正先のGCPの左上 補正先のGCPの右上

補正先のGCPの左下 補正先のGCPの右下

以上のように、補正元の座標値と補正先の座標値をGCPとして取得しておく。

幾何補正の手順その２

OSGeo4W Shellコマンドウィンドウ内で、ogr2ogr -gcpを実行する。
  1. OSGeo4W Shellコマンドウィンドウの起動
    OSGeo4W Shellコマンドウィンドウ
  2. ogr2ogr -gcpの実行
    ファイル名に日本語が含まれた結果
    ファイル名に日本語が含まれるとファイル名として解釈してくれなかった。
    ファイルを開きコマンドが実行されるために、ファイル名を以下のように変更した。
    室蘭駅周辺_tky.shp　→　MuroranEki_tky.shp
    津波_室蘭駅周辺_tky.shp　→　Tsunami_MuroranEki_tky.shp

  3. 再度、ogr2ogr -gcpの実行
    ogr2ogr -gcpを実行
    「Tsunami_MuroranEki_tky.shp」のフィールド名に日本語の「原点X」、「原点Y」が含まれていたため、
    ERROR 1　が出て、「原点X」、「原点Y」が含まれない状態で実行が完了した。
    必要なフィールドが日本語名になっている時は、事前に英数記号に変更が必要。

結果

緑色が、「ogr2ogr -gcp」を実行して幾何補正を行った結果のデータ。
黒色の基図データに合うようにGCPを取得したので、黒色のデータと重なっている。問題無し。
幾何補正の結果

ちなみに「GDAL 1.10.0」が登場するまでは、GRASS 6.4.3 R2（R1では何度試しても途中で強制終了していた）の「幾何補正処理」で幾何補正が可能であった。

ArcGISの「アジャスト」機能に比べると、操作手順など煩わしい。

GRASS 6.4.3 R2と比べると「ogr2ogr -gcp」で操作手順はかなり簡略された。
GRASSは現在、GRASS 6.4.3に更新されているが、操作手順に変更はない模様。

課題

  • 複数ファイルのバッチ処理ができない（GCPの取得を自動処理することが、難題）
  • GCP取得と幾何補正が同一ソフトでできない（ArcGIS、GRASS 6.4.3ならできる）
  • プレビュー確認できない（ArcGISならできる、GRASS 6.4.3はいまいち）

QGIS2.0から「プロセッシング」という表記になった機能がある。
自作したスクリプト処理を実行する事も可能。
「グラフィカルモデラー」で処理をつなげることも可能。
これを利用して上記の課題が解決されることを期待。
大量な処理が出てきたら、社内のプログラマーに依頼し、手っ取り早く作成して貰おうかな。
GDALコマンドでフォルダ内のファイルを幾何補正処理させる事はできるけど。