もうすぐ七夕ですね。みなさんはもう、短冊に「RでもGeoParquetが読みたいです！！！！！」と願い事を書き殴るイメージトレーニングに余念がないことと思います。お疲れ様です。
が、実はもう、RでGeoParquet読めるらしいです。
出典
昨日あった R Consortium のウェビナーの資料を見てて知りました。ここで書くこと以外にも、けっこう知らないこと多そうだったのであとでじっくり見ます。
やり方
やり方は順を追って説明しますが、とりあえず、これが現時点の正解っぽいです。クエリは Overture Maps のドキュメント にあったものを使っています。
library(duckdb)
library(geoarrow)
library(duckplyr)
con <- dbConnect(duckdb::duckdb())
dbExecute(con, "INSTALL spatial")
#> [1] 0
dbExecute(con, "LOAD spatial")
#> [1] 0
dbExecute(con, "INSTALL httpfs")
#> [1] 0
dbExecute(con, "LOAD httpfs")
#> [1] 0
dbExecute(con, "CALL register_geoarrow_extensions()")
#> [1] 0
query <- sql("SELECT
id,
names.primary as name,
class,
geometry -- DuckDB v.1.1.0 will autoload this as a `geometry` type
FROM
read_parquet('s3://overturemaps-us-west-2/release/2025-05-21.0/theme=transportation/type=segment/*', filename=true, hive_partitioning=1)
WHERE
bbox.xmin < 2.314
AND bbox.ymin < 48.882
AND bbox.xmax > 2.276
AND bbox.ymax > 48.865
")
res <- tbl(con, query) |>
arrow::to_arrow() |>
sf::st_as_sf()
res
#> Simple feature collection with 4713 features and 3 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 2.142626 ymin: 48.79561 xmax: 2.453173 ymax: 48.94475
#> CRS: NA
#> First 10 features:
#> id name
#> 1 0841fb47ffffffff046d9d3560f9578c Ligne 18 du Grand Paris Express
#> 2 0841fb47ffffffff047db754a9777dcf Ligne d'Ermont - Eaubonne à Champ-de-Mars
#> 3 08a1fb4629b77fff047ee71c2401d1a7 Villa Jocelyn
#> 4 08b1fb475a4dbfff046def188276a1cd Avenue Victor Hugo
#> 5 08c1fb475a4db3ff046cad1b89c5445f <NA>
#> 6 08c1fb4629b665ff046ee5232c47a412 <NA>
#> 7 08b1fb4629b66fff046cbf214f9de7e3 <NA>
#> 8 08d1fb4629b66cbf046ea71fdb024239 <NA>
#> 9 0891fb475a4fffff046dbf010942869d <NA>
#> 10 08b1fb475a4dbfff046fe51b4b998207 <NA>
#> class geometry
#> 1 <NA> LINESTRING (2.345674 48.917...
#> 2 <NA> LINESTRING (2.287195 48.853...
#> 3 service LINESTRING (2.275723 48.865...
#> 4 secondary LINESTRING (2.275794 48.865...
#> 5 residential LINESTRING (2.276014 48.865...
#> 6 footway LINESTRING (2.276191 48.865...
#> 7 footway LINESTRING (2.276096 48.865...
#> 8 footway LINESTRING (2.276096 48.865...
#> 9 footway LINESTRING (2.275524 48.868...
#> 10 footway LINESTRING (2.275943 48.865...
使うパッケージ
-
duckdb: 今回の主役です。GeoParquet を読みます。 -
sf: R で GIS データと言えばこれ。最終的には sf のデータ形式に変換するのがゴールです。 -
geoarrow: GeoArrow というデータフォーマットを扱うためのパッケージです。DuckDB は GeoArrow 形式でデータを出力できるので、そのデータをこれで読み取るというわけです。 -
duckplyr: DuckDB を dplyr から透過的に扱えるようにするラッパーです。
DuckDB の準備
さて、コードを見ていきましょう。まずはここです。
con <- dbConnect(duckdb::duckdb())
dbExecute(con, "INSTALL spatial")
dbExecute(con, "LOAD spatial")
dbExecute(con, "INSTALL httpfs")
dbExecute(con, "LOAD httpfs")
dbExecute(con, "CALL register_geoarrow_extensions()")
やっていることは、
-
dbConnect()でコネクションを作成する -
spatialとhttpfs（S3 のデータを読むために必要）という拡張を読み込む -
register_geoarrow_extensions()で、データが GeoArrow 形式で吐き出されるようにする
という感じです。
register_geoarrow_extensions() は DuckDB をよく使っている人でも見慣れない関数だと思いますが、これは DuckDB 単体で使うものというより、今回のように DuckDB と他のプログラム（Python や R）が通信するときのために作られたものらしいです。
ちなみに、今はこれを明示的に呼び出さないといけないんですが、デフォルトでこれになってよくない？という提案もあり、将来的には不要になるかもしれません。
クエリを実行
残りは、実際にデータを取得するクエリを実行している部分です。
query <- sql("
...
")
res <- tbl(con, query) |>
arrow::to_arrow() |>
sf::st_as_sf()
この sql() とか tbl() は dbplyr の関数で、まったく dbplyr 的な使い方（dplyr の関数で操作する）はしていないんですが、arrow::to_arrow() が受け付けているデータの形式の関係上？こうなっているみたいです。
たとえば、 DBI には dbGetQueryArrow() という関数があるので「Arrow ってついてるしこれでは？」と思って使ってみたりするんですが、対応していない形式（ちなみに nanoarrow_array_stream というクラスです）らしくエラーになります。
res <- dbGetQueryArrow(con, "...")
res |>
arrow::to_arrow()
#> Error:
#> ! to_arrow() currently only supports Arrow tables, Arrow datasets, Arrow queries, or dbplyr tbls from duckdb connections
また、arrow::to_arrow() の結果は RecordBatchReader なので、こんな感じでいけたりする？とやってみたりもするのですが、これもエラーになります。
arrow::as_record_batch_reader(res) |>
sf::st_as_sf()
#> Error in `st_sf()`:
#> ! no simple features geometry column present
単に RecordBatchReader であるだけではだめで、「どれがジオメトリのカラムなのか」みたいなメタデータが入ってなくては sf::st_as_sf() が変換できない、ということみたいです。ということで、現状、このコードしかなさそうです。
将来の話
いかがだったでしょうか。だいぶスムーズになりましたが、まだちょっと複雑ですよね...。
いま、sf から直接 GeoParquet 読むのは、GDAL 3.11 で入った ADBC ドライバ経由でいけるのでは？という話をしていて、実際わりとできそうな感じなので、たぶん1～2か月後くらいにはもう sf::read_sf() だけでできるようになるような予感がしています。
ということで、続報をお楽しみに！
Discussion