Overture Mapsの建物データ（7億件）からDuckDBのPython APIで日本のデータだけ取得してみた！ #Python

Overture Maps Foundationから、とうとうオープン地図データがリリースされましたね！

世間一般ではどんな反応なのかわかりませんが、界隈では割と盛り上がっている様子です。

Overture Mapsは構造化され、品質が保証された使いやすい地図データをオープンにするべく、整備しているとんでもない規模のプロジェクトで、以下のようなビッグテックが開発を行なっています。

2023/07/26に「Overture 2023-07-26-alpha.0」という名称で初めてデータがリリースされました。
データは4種類に分かれており、個別にライセンスが適用されていますが、基本的には「概ね自由に利用できる」ライセンスになっており、開発者目線だととてもありがたいデータになっています。
（とはいえ、ライセンスはちゃんと確認するようにしましょう。）

場所のテーマ: CDLA Permissive v 2.0
建物のテーマ: ODbL
交通テーマ: ODbL
管理境界のテーマ: CDLA Permissive v 2.0

（https://overturemaps.org/download/より）

色々なデータがありますが、今回は分かりやすい「建物データ」を利用していきましょう。

DuckDBのダウンロード

早速データをダウンロードしていきたいですが、データへのアクセス方法を見ると、どうやらAWSやAzureに保存されているようです。
https://github.com/OvertureMaps/data/blob/main/README.md#how-to-access-overture-maps-data

色々書かれていますが、ちょろっとデータを見るためだけにクラウドを使うのも手間なので、今回は3. DuckDB (SQL)
を利用していきましょう。

DuckDBはOLAP（オンライン上のデータを分散処理で分析すること）に特化したDBであり、OSSです。
軽量かつSQLを利用して分析ができ、Python APIも整備されています。

Pythonによる地理データ分析といえば、GeoPandasなどが挙げられるかと思いますが、より高速に動作するのが特徴のようですね。

詳しくは、以下の記事などが分かりやすかったです。

ということでダウンロードしましょう。

以下でダウンロード可能で、すぐにPythonから利用可能になります。

pip install duckdb

ついでに、その他利用するパッケージも入れておきましょう。

pip install geopandas pyarrow

DuckDBでデータを触ってみる

こちらの手順に従っていきます。

今回はJupyter notebookを利用していきます。
Jupyterの構築方法などは、ググってみてください。

Macユーザーであればこんな記事も書いていますので、参考にしてみてください。

MacでGISデータ分析を始めるためにサクッとAnacondaとjupyter labをインストールしてみる

まずは以下のようにDuckDBのDBを構築します。

import duckdb
import geopandas as gpd
import os

# DB の作成
conn = duckdb.connect("../data/overture.duckdb")
c = conn.cursor()

これで../data/overture.duckdbというファイルが作成されました。

この.duckdb（拡張子は.dbなどでもいいっぽい）自体が、SQLiteのようなデータベース本体になっているようです。

次に、DBにGIS拡張機能をインストールします。。

c.execute(
    """INSTALL spatial;
    INSTALL httpfs;
    LOAD spatial;
    LOAD parquet;
"""
)

データはParquetという、データ分析を高速に行うために特化したファイル形式で提供されています。
（GIS拡張であるGeoParquetについてはこちらで言及しています。
GeoPandasをやるならFlatGeobufより10倍早いGeoParquetを使おう！）

このParquetにはDuckDBから以下のようにアクセスできます。
まずはデータの総数を確認してみましょう。

sql = """
SELECT
    count(*)
FROM
    read_parquet('s3://overturemaps-us-west-2/release/2023-07-26-alpha.0/theme=buildings/type=*/*')
"""

print(c.execute(sql).fetchall())
# [(785524851,)]

7億レコード！！！やっぱり結構多いですね！！

データをダウンロードする

どこに書いていたか忘れてしまいましたが、全データをダウンロードすると150GBを超えるようで、流石にそのまま取り扱うのはめんどくさいので、以下のように、地域を日本のみに絞ってダウンロードしてみましょう。

c.execute(
    """COPY (
    SELECT
           theme,
           type,
           version,
           updateTime,
           JSON(sources) AS sources,
           JSON(names) AS names,
           height,
           numFloors,
           class,
           ST_GeomFromWkb(geometry) AS geometry
      FROM read_parquet('s3://overturemaps-us-west-2/release/2023-07-26-alpha.0/theme=buildings/type=*/*')
     WHERE ST_Within(ST_GeomFromWkb(geometry), ST_Envelope(ST_GeomFromText('POLYGON((122.934570 20.425378, 122.934570 45.551483, 153.986672 45.551483, 153.986672 20.425378, 122.934570 20.425378))')))
) TO '../data/buildings.parquet'
WITH (FORMAT PARQUET);
"""
)

僕の環境では43分かかりましたが、なんとかダウンロードすることができました。

が、これをGeoPandasで読もうとすると…

gdf = gpd.read_parquet("../data/buildings.parquet")
gdf.head()

こんなエラーが出ました。

ValueError: Missing geo metadata in Parquet/Feather file.
            Use pandas.read_parquet/read_feather() instead.

こんなツイートをしている人もいて、どうやらGeoParquetに準拠した形式ではないため、よみこむことができないようです。

仕方ないのでパフォーマンスは落ちますが、FlatGeobufでダウンロードしていきましょう。

c.execute(
    """COPY (
    SELECT
           theme,
           type,
           version,
           updateTime,
           JSON(sources) AS sources,
           JSON(names) AS names,
           height,
           numFloors,
           class,
           ST_GeomFromWkb(geometry) AS geometry
      FROM read_parquet('s3://overturemaps-us-west-2/release/2023-07-26-alpha.0/theme=buildings/type=*/*')
     WHERE ST_Within(ST_GeomFromWkb(geometry), ST_Envelope(ST_GeomFromText('POLYGON((122.934570 20.425378, 122.934570 45.551483, 153.986672 45.551483, 153.986672 20.425378, 122.934570 20.425378))')))
) TO '../data/buildings.fgb'
WITH (FORMAT GDAL, DRIVER 'FlatGeobuf');
"""
)

6時間ほどかかりました…
データ取得方法自体は同じなので、ParquetとFlatGeobufでは書き込み性能に大きな違いがあることが分かります。

ともあれダウンロードができたので、ファイルサイズを確認しましょう。

os.path.getsize("../data/buildings.fgb")
# 7872148968

かなり減りましたが、まだ7GBもあるので、GeoPandasで読み込むにはかなり時間がかかります。
（時間がかかりすぎて今回は諦めました。）

gdf = gpd.read_file("../data/buildings.fgb")
gdf.head()

QGISで読んでみる

今回は大した分析をするつもりはなく、とりあえずデータを取得して閲覧してみよう、くらいの感じですので、ダウンロードしたFlatGeobufはQGISに載せてみましょう！

7GBあるので流石に時間がかかりますが、読み込むと2000万地物まで減っていることが分かります！

だいぶ絞れましたね！

表示するとこんな感じになっていました！

見てみると、東京都でも歯抜けが多いように見えます。

とはいえ、この規模のデータがオープンで公開されるのはとても良いことですね！

おわりに

ということでOverture MapsのデータとDuckDBを使ってデータ取得してみました。

GIS界隈はオープンマインドな方が多く、ツールやデータなど大量に公開されていて、めちゃくちゃ開発がしやすくなっていますし、どんどん便利になってますね。

まだ最初のリリースなので、諸々の問題はあるかもしれませんが、そこはビッグテックの皆様に期待ですね。