🔔 本の下書きプレビューをシェアできるようになりました
Zenn
⛳

e-Stat や国土数値情報の Shapefile を DuckDB で処理するときに詰まったところ（文字コード、CRS）

2025/02/09に公開
GIS
DuckDB
tech

文字コード

e-Stat などで提供されている Shapefile の中には、文字コードが Shift_JIS のものがあります。これを読もうとすると以下のように文字化けになってしまいます。

D SELECT * FROM 'A002005212020XYSWC13/r2ka13.shp' LIMIT 3;
┌──────────┬─────────┬─────────┬─────────┬───────────┬───┬────────────┬───────────┬─────────┬──────────────────────┐
│ KEY_CODE │  PREF   │  CITY   │ S_AREA  │ PREF_NAME │ … │   X_CODE   │  Y_CODE   │ KCODE1  │         geom         │
│ varchar  │ varchar │ varchar │ varchar │  varchar  │   │   double   │  double   │ varchar │       geometry       │
├──────────┼─────────┼─────────┼─────────┼───────────┼───┼────────────┼───────────┼─────────┼──────────────────────┤
│ 13       │ 13      │ 199     │ 001000  │ ?????s         │ … │ 139.822971 │ 35.582316 │ 0010-00 │ POLYGON ((192.3304…  │
│ 13       │ 13      │ 402     │ 000000  │ ?????s         │ … │ 139.917477 │ 31.887491 │ 0000-00 │ POLYGON ((8173.880…  │
│ 13       │ 13      │ 402     │ 000000  │ ?????s         │ … │ 140.051022 │ 31.438999 │ 0000-00 │ POLYGON ((20721.04…  │
├──────────┴─────────┴─────────┴─────────┴───────────┴───┴────────────┴───────────┴─────────┴──────────────────────┤
│ 3 rows                                                                                      30 columns (9 shown) │
└──────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘

DuckDB は CSV の読み込みでも UTF-8 にしか対応していなかった[1]ので、方法なさそうかな...と思って天を仰いだのですが、よく考えてみると、DuckDB が内部的に使っている GDAL は文字コードの指定ができるので、DuckDB にもなんかオプションあるのでは...？と思って調べたら、ありました。

ST_Read() に open_options = ['ENCODING=CP932'] を指定するといいみたいです。

D SELECT * FROM ST_Read('A002005212020XYSWC13/r2ka13.shp', open_options = ['ENCODING=CP932']) LIMIT 3;
┌──────────┬─────────┬─────────┬─────────┬───────────┬───┬────────────┬───────────┬─────────┬──────────────────────┐
│ KEY_CODE │  PREF   │  CITY   │ S_AREA  │ PREF_NAME │ … │   X_CODE   │  Y_CODE   │ KCODE1  │         geom         │
│ varchar  │ varchar │ varchar │ varchar │  varchar  │   │   double   │  double   │ varchar │       geometry       │
├──────────┼─────────┼─────────┼─────────┼───────────┼───┼────────────┼───────────┼─────────┼──────────────────────┤
│ 13       │ 13      │ 199     │ 001000  │ 東京都    │ … │ 139.822971 │ 35.582316 │ 0010-00 │ POLYGON ((192.3304…  │
│ 13       │ 13      │ 402     │ 000000  │ 東京都    │ … │ 139.917477 │ 31.887491 │ 0000-00 │ POLYGON ((8173.880…  │
│ 13       │ 13      │ 402     │ 000000  │ 東京都    │ … │ 140.051022 │ 31.438999 │ 0000-00 │ POLYGON ((20721.04…  │
├──────────┴─────────┴─────────┴─────────┴───────────┴───┴────────────┴───────────┴─────────┴──────────────────────┤
│ 3 rows                                                                                      30 columns (9 shown) │
└──────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘

ドキュメントはこのあたりです。

https://duckdb.org/docs/extensions/spatial/functions.html#st_read

CRS

DuckDB の spatial extension は ST_Transform() などで自在に投影変換できますが、DuckDB 上のデータには CRS が保持されていません。なので、自分で CRS を指定する必要があるみたいです。具体的には、

  • ST_Transform() の source CRS
  • COPY TO で書き出すときの SRS オプション

というあたりです。

ST_Transform()

比較のために、まずは PostGIS を見てみましょう。PostGIS の ST_Transform() のシグネチャ（ドキュメント）を見ると、以下のように、from_proj があるバージョンとないバージョンがあるのがわかると思います。from_proj を明示的に指定しない場合は、そのカラムの CRS が source として使われます。

geometry ST_Transform(geometry g1, integer srid);
geometry ST_Transform(geometry geom, text to_proj);
geometry ST_Transform(geometry geom, text from_proj, text to_proj);
geometry ST_Transform(geometry geom, text from_proj, integer to_srid);

一方、DuckDB の ST_Transform() のシグネチャ（ドキュメント）を見ると、すべて source_crs が入っているのがわかると思います。つまり、ST_Transform() を使うには、ユーザーがそのデータの CRS をあらかじめ知っておく必要があります。

BOX_2D ST_Transform (geom BOX_2D, source_crs VARCHAR, target_crs VARCHAR)
BOX_2D ST_Transform (geom BOX_2D, source_crs VARCHAR, target_crs VARCHAR, always_xy BOOLEAN)
POINT_2D ST_Transform (geom POINT_2D, source_crs VARCHAR, target_crs VARCHAR)
POINT_2D ST_Transform (geom POINT_2D, source_crs VARCHAR, target_crs VARCHAR, always_xy BOOLEAN)
GEOMETRY ST_Transform (geom GEOMETRY, source_crs VARCHAR, target_crs VARCHAR)
GEOMETRY ST_Transform (geom GEOMETRY, source_crs VARCHAR, target_crs VARCHAR, always_xy BOOLEAN)

DuckDB 側でも、このために ST_Read_Meta() という関数が用意されました。ちょっとひと手間という感じではありますが、これを使えば DuckDB の中で作業を完結できます。

https://github.com/duckdb/duckdb-spatial/issues/156

具体的には、以下のようなクエリを実行すると CRS が見れます（複数の CRS が混在している場合もあるので array の中に入っていますが、ここでは、CRS はひとつしかないと考えて [1] 決め打ちで取り出しています）。このデータでいうと、'EPSG:2451' を指定することになります。

D SELECT unnest(layers[1].geometry_fields[1].crs) FROM St_read_meta('A002005212020XYSWC13/r2ka13.shp');
┌──────────────────────┬───────────┬───────────┬──────────────────────┬──────────────────────┬──────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│         name         │ auth_name │ auth_code │         wkt          │        proj4         │                                                     projjson                                                     │
│       varchar        │  varchar  │  varchar  │       varchar        │       varchar        │                                                     varchar                                                      │
├──────────────────────┼───────────┼───────────┼──────────────────────┼──────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ JGD2000 / Japan Pl…  │ EPSG      │ 2451      │ PROJCS["JGD2000 / …  │ +proj=tmerc +lat_0…  │ {\n  "$schema": "https://proj.org/schemas/v0.5/projjson.schema.json",\n  "type": "ProjectedCRS",\n  "name": "J…  │
└──────────────────────┴───────────┴───────────┴──────────────────────┴──────────────────────┴──────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘

COPY TO

COPY ... TO 'path/to/file' WITH (FORMAT GDAL, ...) で書き出す際にも、CRS を指定する必要があります。例えば、Shapefile で言うと、読み込んだものをそのまま書き出して元のファイルと比べると、.prj ファイルがなくなっていることに気付くと思います。

CREATE TABLE poly AS SELECT * FROM ST_Read('A002005212020XYSWC13\r2ka13.shp', open_options = ['ENCODING=CP932']);
COPY poly TO 'foo.shp' WITH (FORMAT GDAL, DRIVER 'ESRI Shapefile');

結果：

ではどうすればいいかというと、SRS というオプションを指定します（ドキュメント）。こんな感じです。

COPY poly TO 'foo2.shp' WITH (FORMAT GDAL, DRIVER 'ESRI Shapefile', SRS 'EPSG:2451');

ここでは .prj ファイルの有無が見た目で分かりやすかったので Shapefile にしましたが、
もちろん Shapefile 以外のフォーマットでも同じように指定できます。例えば、Web メルカトルに投影変換してから GeoPackage 形式で保存するには、以下のような SQL を実行します。

UPDATE poly SET geom = ST_Transform(geom, 'EPSG:2451', 'EPSG:3857');

COPY poly TO 'foo.gpkg' WITH (FORMAT GDAL, DRIVER 'GPKG', SRS 'EPSG:3857');

今後の計画？

詳細はちょっと今の私では理解できなかったのですが、カラムの型情報に CRS を持たせる？という計画はあるようです。これができれば、自分で元の CRS を意識する必要はなくなるのかもしれません。

https://github.com/duckdb/duckdb-spatial/issues/441

脚注

  1. v1.2 で Latin-1 と UTF-16 はサポートされました ↩︎

Discussion

ログインするとコメントできます