e-Stat や国土数値情報の Shapefile を DuckDB で処理するときに詰まったところ(文字コード、CRS)
文字コード
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) │
└──────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
ドキュメントはこのあたりです。
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 の中で作業を完結できます。
具体的には、以下のようなクエリを実行すると 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 を意識する必要はなくなるのかもしれません。
-
v1.2 で Latin-1 と UTF-16 はサポートされました ↩︎
Discussion