Mooの技術メモ

技術的なメモ中心。

Overture MapsのデータをRで可視化する

R GIS 地図

Overture Mapsのデータが先月公開されたので軽く可視化してみる。 overturemaps.org

Parquet形式で提供されているようで、以下のツールを使ってデータを取得するサンプルがGitHubで公開されている。

  1. Amazon Athena
  2. Microsoft Synapse
  3. DuckDB

3つめのDuckDBはRのパッケージがあるようだ。

最近Rを触っていないしとりあえずローカルの環境で完結するので今回はDuckDBを使う。

使用したパッケージ

DuckDBに接続し、データをダウンロードする

データベースのインスタンスを作成し接続する。

con = dbConnect(duckdb(), dbdir = ":memory:")

duckdb() でDBのインスタンスを作成するが、dbdir=":memory:" とすることでメモリ上に作成することができる。 今回はデータをダウンロードしたらそれ以降は使わないので、メモリ上に作成することにした。

次のコードでデータをローカルにコピーするSQLを実行する。

dbExecute(
    con,
    "
    INSTALL spatial;
    INSTALL httpfs;
    INSTALL json;
    LOAD spatial;
    LOAD httpfs;
    LOAD json;
    COPY (
        SELECT
            type,
            subType,
            localityType,
            adminLevel,
            isoCountryCodeAlpha2,
            JSON(names) AS names,
            JSON(sources) AS sources,
            ST_GeomFromWkb(geometry) AS geometry
            FROM read_parquet('s3://overturemaps-us-west-2/release/2023-07-26-alpha.0/theme=admins/type=*/*', filename=true, hive_partitioning=1)
            WHERE adminLevel = 2
            AND ST_GeometryType(ST_GeomFromWkb(geometry)) IN ('POLYGON','MULTIPOLYGON')
        ) TO 'countries.json'
    WITH (FORMAT GDAL DRIVER 'GeoJSON');
    "
)

dbDisconnect(con)

DuckDBのエクステンションが必要になるので、spatial、httpfs、jsonの3つをインストール・ロードしている。 クエリはGitHubで公開されているサンプルをそのまま使うことにした。

サンプルではフォーマットをParquetではなくGeoJSONにしているので、これParquetでダウンロードしたらGeoParquetじゃないってオチだろ、と思ったら案の定そのとおりでGitHubにdiscussionが上がっていた。

github.com

そのうちGeoParquetでDLできるようになるとは思うが、今後に期待。

ダウンロードができたらもう不要なのでDBとの接続は切っておく。

GeoJSONを読み込んで可視化

次のコードでデータをsf形式で読み込む。 面積だけ計算し、使わなさそうなカラムは削除した。

countries <- st_read("countries.json") %>% 
    mutate(area = st_area(.$geometry)) %>% 
    mutate(area_km2 = set_units(area, km2)) %>% 
    select(isocountrycodealpha2, area_km2)

可視化にはtmapを使って面積でchoropleth mapを作った。

tmap_mode("view")
tm_shape(countries) +
    tm_fill("area", palette = "Blues", alpha = 0.8) +
    tm_borders(col = "grey40", lwd = 0.4, lty = "solid") +
    tm_layout(title = "Overture Maps countries data", legend.title.size = 1, legend.text.size = 0.6) +
    tm_view(view.legend.position = c("left", "bottom"))

カナダが雑すぎるw

国別に塗り分けたマップも作成してみた。

tm_shape(countries) +
    tm_polygons("isocountrycodealpha2",
                alpha = 0.8,
                title = "Country code",
                popup.vars = c("Country code" = "isocountrycodealpha2")) +
    tm_style("classic") +
    tm_layout(title = "Overture Maps countries data", legend.title.size = 1, legend.text.size = 0.6) +
    tm_view(view.legend.position = c("right", "bottom"))

四国と本州が癒着してる。氷河期かな。

Country codeはISO 3166-1 alpha-2。北方領土みたいな係争地はuser assigned codeになってるっぽい。