はじめに
「GISエンジニア」というデカい主語で始まった本記事。最近話題のDuckDBですが、GISデータを扱うのに非常に有用で、今までGeoPandasやPostGISでやっていたような処理を、CLIで・最小の環境構築で・高速に実行できます。本記事では、GISデータのプロセッシングという観点で、使い方を紹介します。
DuckDBとは何か
正確な定義については私は間違えそうなので公式Webサイトなどをご覧いただくとして。
「DBというくらいだから、データベースなんでしょ？」正解です。しかしPostgreSQLなどのようなある種「重厚な」DBソフトウェアとはちょっと違います。PostgreSQLでいうとpsqlクライアントだけが利用出来るようなイメージです（厳密には永続的なデータベースも作成は出来るが）。
私の理解では、DuckDBは「DBクライアント」であり「CLIツール」です。DuckDBクライアントから、様々なデータベースに接続し・SQLで操作できる。これに尽きます。
DuckDBがユニークなのは、PostgreSQLやMySQLといった伝統的なRDBを操作出来ることではなく、たとえばCSVファイルやXLSXを開いてSQLを実行出来たり、リモートファイルにHTTPS越しで接続したり、あるいはAWS S3などのオブジェクトストレージに対してスキャンのようなことが出来たり、「実にさまざまなデータを開いてSQLで操作できる」という点です。パフォーマンスもよく軽快に動作します。
Extension
DuckDBはExtensionという機構を備えており、後付けで新しいデータ形式をサポートしたり出来ます。サードパーティ向けの機能でありつつも、重要度の高いExtensionは公式ドキュメントに列挙されています。
https://duckdb.org/docs/extensions/overview
たとえばHTTP越しのデータ操作についてはhttpfsなんていうExtensionを利用します。本記事ではspatialを利用してGISデータを操作する話をします。
DuckDBのインストール
https://duckdb.org/docs/installation/
公式Webサイトで実行環境別にインストール方法を示してくれています。Windowsでもネイティブ動作するのもうれしい人が多いかもしれない。インストールがすぐ完了するのもいいところですね。
spatial Extension
https://duckdb.org/docs/extensions/spatial/overview
このExtensionは以下の機能を提供します：
- GEOMETRY型のサポート・各種関数群の提供
- 任意形式のGISベクターデータの読み書き（ラスターは非対応）
- OSMデータの読み込み（
osm.pbfを直接、表形式で読み込める） - 空間インデックスの構築（比較的最近実装された）
かみ砕くと、「GISデータを手軽に表形式で開くことが出来るようになり」「SQLでデータを整形でき」「必要に応じて空間関数を利用し」「任意形式で出力できる」ということです。
SQLでデータを処理したいからと言ってPostGISにデータを投入して…というのはつらいものがあります。多分そういうケースではGeoPandasやQGISあたりで処理することが多いのかなと思いますが、DuckDBならコマンドラインですぐに立ち上げて、GISデータを直接ロード・編集・書き出しできます。
DuckDB + Spatialの使い方
DuckDBを起動します。
duckdb
# これでオンメモリで起動する。2つ目の引数でファイル名を指定するとデータベースファイルが作成される。
spatial Extensionを利用するには以下のSQLを実行する必要があります。
INSTALL spatial; -- 一度実行すれば、次回以降は不要
LOAD spatial; -- クライアントを起動するたびに、実行する必要がある
データは、例として「国土数値情報 - 行政区域データ」を利用します。
https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-2024.html
GISデータを読み込む
なんとこんなSQLでデータを表形式に展開できます
SELECT * FROM 'N03-20240101_GML/N03-20240101.shp';
┌─────────┬────────────────┬──────────┬──────────┬─────────┬─────────┬─────────────────────────────────────────────────┐
│ N03_001 │ N03_002 │ N03_003 │ N03_004 │ N03_005 │ N03_007 │ geom │
│ varchar │ varchar │ varchar │ varchar │ varchar │ varchar │ geometry │
├─────────┼────────────────┼──────────┼──────────┼─────────┼─────────┼─────────────────────────────────────────────────┤
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 中央区 │ 01101 │ POLYGON ((141.25694 42.997815279, 141.2563044… │
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 北区 │ 01102 │ POLYGON ((141.333254449 43.075053892, 141.332… │
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 東区 │ 01103 │ POLYGON ((141.373410558 43.068399721, 141.373… │
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 白石区 │ 01104 │ POLYGON ((141.38201917 43.048316387, 141.3819… │
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 豊平区 │ 01105 │ POLYGON ((141.363712866 42.941538225, 141.363… │
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 南区 │ 01106 │ POLYGON ((141.235418054 42.822725279, 141.235… │
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 西区 │ 01107 │ POLYGON ((141.169036109 43.082901387, 141.170… │
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 厚別区 │ 01108 │ POLYGON ((141.504765551 43.021098892, 141.504… │
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 手稲区 │ 01109 │ POLYGON ((141.227041946 43.083017505, 141.226… │
│ 北海道 │ 石狩振興局 │ │ 札幌市 │ 清田区 │ 01110 │ POLYGON ((141.413818612 42.89511, 141.4140658… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((140.719181582 41.749965748, 140.719… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((140.711850726 41.770813505, 140.714… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((141.030827782 41.716509441, 141.030… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((140.701418521 41.787353586, 140.701… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((140.69675703 41.792666586, 140.6929… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((140.949656109 41.937226667, 140.949… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((141.028324163 41.717673054, 141.028… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((141.007578612 41.714451387, 141.007… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((141.185812361 41.804415081, 141.185… │
│ 北海道 │ 渡島総合振興局 │ │ 函館市 │ │ 01202 │ POLYGON ((140.818869805 41.764438694, 140.818… │
│ ・ │ ・ │ ・ │ ・ │ ・ │ ・ │ ・
│
│ ・ │ ・ │ ・ │ ・ │ ・ │ ・ │ ・
│
│ ・ │ ・ │ ・ │ ・ │ ・ │ ・ │ ・
│
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((123.042224514 24.460154964, 123.042… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.998083489 24.439785946, 122.998… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((123.042094488 24.460186099, 123.042… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.998920506 24.439903072, 122.998… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.998780674 24.439765459, 122.998… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.997914981 24.439807216, 122.997… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.942213074 24.443955297, 122.942… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.994233476 24.440338189, 122.994… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((123.038420571 24.458696685, 123.038… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((123.042255577 24.460196243, 123.042… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.998532322 24.439807423, 122.998… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.99277716 24.440740721, 122.9927… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.995011907 24.440169126, 122.995… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((123.005067445 24.439328811, 123.005… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.934714436 24.447848991, 122.934… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.932986602 24.451701342, 122.932… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.933963606 24.44820464, 122.9339… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.942131634 24.443952748, 122.942… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.998760636 24.439854613, 122.998… │
│ 沖縄県 │ │ 八重山郡 │ 与那国町 │ │ 47382 │ POLYGON ((122.933082724 24.451764126, 122.933… │
├─────────┴────────────────┴──────────┴──────────┴─────────┴─────────┴─────────────────────────────────────────────────┤
│ 124133 rows (40 shown) 7 columns │
└──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
シェープファイルが持つ各カラムのほかに、geomカラムがGEOMETRY型として含まれていることに留意してください。shp,gpkb,fgbは拡張子から自動的にspatialで読み込まれます。それ以外の場合はST_Readを利用します。行政区域データにはGeoJSONファイルが同梱されているので：
SELECT * FROM ST_Read('N03-20240101_GML/N03-20240101.geojson');
これを実行すると、先ほどと同じような結果が得られます。なおSHP形式よりは読み込みが遅いであろうことに注意を払いましょう。フォーマットごとのメリット・デメリットは以下の記事を参照してください。
https://qiita.com/Kanahiro/items/1cc64aecf5cf9d5dd2a0#geojson
GISデータの処理
慣れ親しんだSQLの文法がおおむね利用できるはずなので、思い思いのSQLで君だけのGISデータを作ろう！
CREATE TABLE admin AS SELECT * FROM 'N03-20240101_GML/N03-20240101.shp';
-- 北海道だけを抽出してみる
CREATE TABLE hokkaido AS SELECT * FROM admin WHERE N03_001 = '北海道';
-- 市町村を面積の大きい順に並べてみる
SELECT N03_004, SUM(ST_Area(geom)) as area FROM hokkaido GROUP BY N03_004 ORDER BY area
DESC;
┌────────────┬──────────────────────┐
│ N03_004 │ area │
│ varchar │ double │
├────────────┼──────────────────────┤
│ 留別村 │ 0.16813080072817768 │
│ 北見市 │ 0.1597763516964302 │
│ 足寄町 │ 0.15638075224608278 │
│ 釧路市 │ 0.15087898282048684 │
│ 別海町 │ 0.1493826026017987 │
│ 遠軽町 │ 0.14920640301425392 │
│ 枝幸町 │ 0.12700554764561758 │
│ 士別市 │ 0.12582808735599668 │
│ 新ひだか町 │ 0.12553869682247304 │
│ 札幌市 │ 0.12376820858364082 │
│ 標茶町 │ 0.12193811855169379 │
│ 新得町 │ 0.11804222416382357 │
│ 上川町 │ 0.11727692183000576 │
│ 日高町 │ 0.10904120425577404 │
経緯度平面での計算となっています。 球面（PostGISでいうGEOGRAPHY）を考慮するには以下のとおりです。
SELECT N03_004, SUM(ST_Area_Spheroid(ST_FlipCoordinates(geom))) as area FROM hokkaido GROUP BY N03_004 ORDER BY area DESC;
┌────────────┬────────────────────┐
│ N03_004 │ area │
│ varchar │ double │
├────────────┼────────────────────┤
│ 留別村 │ 1477507333.2948418 │
│ 北見市 │ 1427412240.979009 │
│ 足寄町 │ 1408042799.9170861 │
│ 釧路市 │ 1363255071.8516328 │
│ 別海町 │ 1344657747.2254722 │
│ 遠軽町 │ 1331769686.1099417 │
│ 新ひだか町 │ 1147553831.5218844 │
│ 札幌市 │ 1121261297.1618342 │
│ 士別市 │ 1119215025.4770842 │
│ 枝幸町 │ 1117183402.6681747 │
│ 標茶町 │ 1099369005.6460247 │
│ 新得町 │ 1063832711.8133336 │
│ 上川町 │ 1049472809.7730615 │
│ 日高町 │ 992065308.4878501 │
新たにジオメトリを生成することもできます。ポリゴンの内部保証点のテーブルを作ってみましょう。
CREATE TABLE poi AS SELECT ST_PointOnSurface(geom) as geom, N03_004 as name, N03_001 as pref FROM hokkaido;
GISデータの書き出し
加工したテーブルをGISデータのファイルとして出力してみましょう。
COPY poi TO 'poi.gpkg' WITH (FORMAT GDAL, DRIVER 'GPKG');
このSQLだけで、poiテーブルがGeoPackage形式で出力されます。QGISで確認してみると：
正常にポイントが生成されていることがわかりますね
COPYはSELECTと組み合わせても利用できるので
COPY (SELECT * EXCLUDE(geom), ST_Simplify(geom, 0.1) FROM hokkaido) TO 'simplified.gpkg' WITH (FORMAT GDAL, DRIVER 'GPKG');
過度に抽象化されてしまったhokkaidoくん
:::note info QGISというと、DuckDBデータベースからジオメトリを直接読み出せるQDuckDBというプラグインも公開されています。
https://oslandia.gitlab.io/qgis/qduckdb/ :::
終わりに
DuckDBを用いたGISデータの読み込み・加工・書き出しの一連の流れを紹介しました。いかがでしたでしょうか。明日から使いたくなりそうな便利なソフトウェアですよね。GISに限らないいろいろな使い道があるので、ぜひまずはGISデータの編集ツールとして使い始めてみてください。すぐに深堀したくなると思います。