久々の投稿です。GISツール等で使うポリゴンのデータ加工をするため、GDALのogr2ogr ツールをどう使えるかを解説した記事となります。
qgis-japan-mesh プラグイン
わたしが所属する某M社ではOSSのGISツール・QGIS向けのオープンソースのプラグインをいくつか開発しており、その中に qgis-japan-mesh というプラグインがあります。
これは地域メッシュコードと呼ばれるメッシュをQGISで扱いたい場合、便利なプラグインです。
qgis-japan-mesh プラグインの使い方をざっと説明します。
プラグインを起動し、お好きな地域メッシュ区画を選択して[実行]。
その後、しばらく待ちます。
処理が完了すると、日本全体に網のようなものがかかります。
これをズームしてみると、
細かいメッシュ状のポリゴンで囲まれていることが分かります。
これが地域メッシュです。
主に統計などに使う目的で、緯度・経度に基づき日本の地域を隙間なく網の目の広さの区画に分けたもので、扱う範囲により、1次～3次、さらに細かなものがあります。
メッシュごとに地域メッシュコードと呼ばれる番号が設定されており、番号からその区画の緯度経度もおおよそ分かります。
QGISでは地域メッシュを都道府県に合わせて「美しく」切り取れる？
このように qgis-japan-mesh プラグインを使うことで、簡単に日本全体を囲む地域メッシュを作ることはできました。
ただ、日本全体の地域メッシュではなく、一部の都道府県の範囲に含まれる地域メッシュが欲しかったとします。
ということで、サンプルで扱う秋田県のポリゴンを用意します。
結果が分かりやすいため、地域メッシュは格子間隔が10kmと広めの２次メッシュを使います。
ここから、秋田県のポリゴンの形に含まれるメッシュ群が欲しいですね。
QGISのメニューから [ベクタ] > [空間演算ツール] > [交差（intersect）] を選択してみます。
- 入力レイヤ: 2次メッシュのレイヤー
- オーバーレイレイヤ: 都道府県のポリゴン表示レイヤー
これで [実行] してみると、こんな絵が出てきます。
さて、このようにして出てきた絵をみて、わたしは思いました。
「欲しかった絵はコレではなかった。美しくない」
QGIS空間演算ツールの [交差（intersect）] では、2次メッシュポリゴンのレイヤーにオーバーレイして重ねた都道府県のポリゴンの外枠に合わせ、メッシュの格子自体の形も削ってしまいます。
でも、わたしが欲しかったのは、都道府県のポリゴンを囲んだメッシュ格子の集合。
文章だと分かりづらいので、ざっくりとした絵で描いて表現してみると、こんな感じのカクカクした結果が欲しかったのです。
さて、こういったメッシュを得るためには、地域メッシュに対し、どのような加工を行えばいいでしょうか？
ogr2ogr を使って、地域メッシュを都道府県に合わせて「美しく」切り取ってみる
ということで、データの加工に GDAL ツールに含まれる ogr2ogr を使ってみたいと思います。
まず、前準備として、加工するレイヤー同士をファイルに保存が必要です。
QGISで表示しているレイヤーから、GPKGファイルにエクスポートします。
- ２次メッシュのレイヤーをGPKGにエクスポート ... (1) akitaken.gpkg
- 都道府県のポリゴンのレイヤーもGPKGにエクスポート ... (2) 2nd_mesh.gpkg
次に、GDALのogr2ogrを使い、エクスポートした二つのファイルに対し、以下のコマンドを実行します。
ogr2ogr -update -append 2nd_mesh.gpkg akitaken.gpkg -nln akitaken
ogr2ogr -f GPKG 2nd_mesh_cliped.gpkg -dialect SQLITE -sql "SELECT A.* FROM '2' as A, 'akitaken' as B WHERE ST_Intersects(A.geom, B.geom)" 2nd_mesh.gpkg
この２つのコマンドを実行することで 2nd_mesh_cliped.gpkg ファイルが生成されました。
これを QGIS で開いてみるとこのような絵が出てきます。
どうやらこれで、わたしがイメージしてた絵を作ることができました。
やったね。
さて、２つのGPKGファイルを作り、２つのコマンドを実行することで、あっさりとわたしがイメージ通りのポリゴンを取得することができましたが、ここまででどういった加工を行ったのか、詳しく解説してみます。
実行したコマンドの意味
GPKGファイルについて
GPKGファイルは、Geopackageというファイル形式です。
QGIS ver 3.0以降でサポートするようになったGIS用のファイル形式で、データ構造としてはポリゴン・線・点等の位置情報で用いるジオメトリ列をサポートするSQLiteのデータベースです。
したがって、SQLiteに対しSQLクエリを発行できるツールを用いてデータ抽出することも可能です。
OGRツールについて
ポリゴン・線・点はGIS的にはベクトルデータです。GDALでベクトルデータを扱うにはOGRツールを用います。
QGISのレイヤーをエクスポートしたファイルについて、ogrinfo コマンドでファイルのメタ情報を覗いてみます。
$ ogrinfo akitaken.gpkg
INFO: Open of `akitaken.gpkg'
using driver `GPKG' successful.
1: todoufukenn__output (Multi Polygon)
$ ogrinfo 2nd_mesh.gpkg
INFO: Open of `2nd_mesh.gpkg'
using driver `GPKG' successful.
1: 2 (Polygon)
akitaken.gpkgは、todoufukenn__output (Multi Polygon) 2nd_mesh.gpkgは、2 (Polygon) と表示されました。
これは、それぞれのファイルでジオメトリ列を含むテーブルの名前と、ジオメトリ列のデータ種類を示します。このテーブルは、QGISではレイヤーとして用いることができます。
最初の ogr2ogr コマンドでは何をやってる？ > 加工するファイルのレイヤーにオーバーレイさせたいレイヤーを追加している
次に実行した ogr2ogr のコマンドを見てみます。最初に発行したコマンドはこのようなものでした。
ogr2ogr -update -append 2nd_mesh.gpkg akitaken.gpkg -nln akitaken
このコマンドは、akitaken.gpkgが含んでいるテーブル todoufukenn__output を、2nd_mesh.gpkg に追加しています。ついでに、-nln オプションによって、テーブル名を akitaken と設定しています。
今回は一つのSQLiteデータベースに含んでいる2つのテーブルを結合した結果を抽出したい、すなわち、データの加工元レイヤーにオーバーレイするレイヤーを追加したく、このような処理を行っています。
2番目の ogr2ogr コマンドでは何をやっているか？ > レイヤー同士をSQLでJOIN。条件に一致する結果を抽出し、結果を別ファイルに保存する
次、2番目に実行した ogr2ogr コマンドでは何をやっているのでしょうか？
ogr2ogr -f GPKG 2nd_mesh_cliped.gpkg -dialect SQLITE -sql "SELECT A.* FROM '2' as A, 'akitaken' as B WHERE ST_Intersects(A.geom, B.geom)" 2nd_mesh.gpkg
ここでは、2つのテーブルを含む 2nd_mesh.gpkg に対し、 -dialect SQLITE -sql [SQL文] というオプションにより、指定したSQLで2nd_mesh.gpkgに含まれるテーブルからデータを抽出しています。
そして、抽出した結果を -f GPKG 2nd_mesh_cliped.gpkg というオプションにより、2nd_mesh_cliped.gpkg というGeoPackageファイルに保存します。
SQLだけ取り出すとこのようなものです。
SELECT A.*
FROM '2' as A, 'akitaken' as B
WHERE ST_Intersects(A.geom, B.geom)
2nd_mesh.gpkg に含んでいる '2'テーブル と 'akitaken' テーブル を結合。
結合条件はWHERE文の ST_Intersects(A.geom, B.geom) となります。
ST_Intersects 関数は、ogr2ogrに-dialect SQLITEを指定した時に使える空間演算関数ですが、標準のSQLiteがサポートする関数ではありません。
Geopackageはジオメトリを扱うことができるSQLite互換のデータベースと説明しましたが、こちらでサポートしている空間演算関数となります。
関数の機能は、2つの引数で指定した各テーブルのジオメトリ同士が交差しているレコードだけを抽出するというものです。
今回は、全国分の2次メッシュテーブルのジオメトリから、秋田県のポリゴンに含まれる結果だけを抽出しています。
そしてこの条件を設定したSELECT文による抽出列は、2次メッシュを管理してる"2"のテーブル結果だけです。
これはつまり、ふたつのポリゴン同士が重なる結果から、秋田県のポリゴンに含まれる、未加工の2次メッシュ一覧だけを抽出してるということです。カクカクした結果が得られるのはこういう訳なのですね。
おわりに
わたしはQGISよりはSQLの方が得意なので、このようなやり方になりました。
もしかすると、QGISだけでも今回の結果は抽出することができるかもしれませんが、あいにくとやり方が分かりません。
もしご存じの方がいれば、コメント等でお知らせいただけると幸いです。
やり方をこうすればいいのでは？というアイデアは同僚の @KEI_YAMA さんからいただいたものを使ってます。ありがとうございました。