GISエンジニアの必携ツールとなったDuckDBの紹介

# はじめに

「GISエンジニア」というデカい主語で始まった本記事。最近話題のDuckDBですが、GISデータを扱うのに非常に有用で、今までGeoPandasやPostGISでやっていたような処理を、CLIで・最小の環境構築で・高速に実行できます。本記事では、GISデータのプロセッシングという観点で、使い方を紹介します。

# DuckDBとは何か

https://duckdb.org/

正確な定義については私は間違えそうなので公式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を起動します。

```sh
duckdb
# これでオンメモリで起動する。2つ目の引数でファイル名を指定するとデータベースファイルが作成される。
```

`spatial` Extensionを利用するには以下のSQLを実行する必要があります。

```sql
INSTALL spatial; -- 一度実行すれば、次回以降は不要
LOAD spatial; -- クライアントを起動するたびに、実行する必要がある
```

データは、例として「国土数値情報 - 行政区域データ」を利用します。

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-2024.html

## GISデータを読み込む

なんとこんなSQLでデータを表形式に展開できます

```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ファイルが同梱されているので：

```sql
SELECT * FROM ST_Read('N03-20240101_GML/N03-20240101.geojson');
```

これを実行すると、先ほどと同じような結果が得られます。なおSHP形式よりは読み込みが遅いであろうことに注意を払いましょう。フォーマットごとのメリット・デメリットは以下の記事を参照してください。

https://qiita.com/Kanahiro/items/1cc64aecf5cf9d5dd2a0#geojson

## GISデータの処理

慣れ親しんだSQLの文法がおおむね利用できるはずなので、思い思いのSQLで君だけのGISデータを作ろう！

```sql
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）を考慮するには以下のとおりです。

```sql
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 │
```

新たにジオメトリを生成することもできます。ポリゴンの内部保証点のテーブルを作ってみましょう。

```sql
CREATE TABLE poi AS SELECT ST_PointOnSurface(geom) as geom, N03_004 as name, N03_001 as pref FROM hokkaido;
```

## GISデータの書き出し

加工したテーブルをGISデータのファイルとして出力してみましょう。

```sql
COPY poi TO 'poi.gpkg' WITH (FORMAT GDAL, DRIVER 'GPKG');
```

このSQLだけで、`poi`テーブルがGeoPackage形式で出力されます。QGISで確認してみると：

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/73197/1aa38eb1-5523-c7a0-e442-9feb0f05d6b0.png)
*正常にポイントが生成されていることがわかりますね*

`COPY`は`SELECT`と組み合わせても利用できるので

```sql
COPY (SELECT * EXCLUDE(geom), ST_Simplify(geom, 0.1) FROM hokkaido) TO 'simplified.gpkg' WITH (FORMAT GDAL, DRIVER 'GPKG');
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/73197/c32686d8-8297-8657-0a9f-429dbd45ccf4.png)
*過度に抽象化されてしまったhokkaidoくん*


:::note info
QGISというと、DuckDBデータベースからジオメトリを直接読み出せるQDuckDBというプラグインも公開されています。

https://oslandia.gitlab.io/qgis/qduckdb/
:::

## 終わりに

DuckDBを用いたGISデータの読み込み・加工・書き出しの一連の流れを紹介しました。いかがでしたでしょうか。明日から使いたくなりそうな便利なソフトウェアですよね。GISに限らないいろいろな使い道があるので、ぜひまずはGISデータの編集ツールとして使い始めてみてください。すぐに深堀したくなると思います。