地球は平らではなく、あなたのボロノイ図もそうであるべきではありません

美とファッションに関する地球は平らではなく、あなたのボロノイ図も同じようにすべきではありません

精度に関する物語、Pythonで球体ジオスパイアル・ボロノイ図の力を明らかにする

あなたはおそらくボロノイ図とその地球空間解析での使用について聞いたことがあるかもしれません。もし聞いたことがなければ、ここで簡単に説明します：ボロノイ図は、与えられたシードよりも他のどの点よりも平面上のすべての点が近い領域に平面を分割します。それは数学者ゲオルギー・ボロノイにちなんで名付けられています。詳しくはWikipediaをご覧ください。

これが地球空間ドメインにどのように適用されるのでしょうか？ボロノイ図を使用することで、大きなスケールでのある都市の住民にとって最寄りの公共交通停留所を個々の建物ごとに計算するよりも、より速く見つけることができます。また、異なるブランド間の市場シェア分析にも使用することができます。

この記事では、平面上の投影された座標で計算された通常のボロノイ図と球体のボロノイ図との違いを示し、後者の優位性を示したいと思います。

寸法と投影 – なぜ重要なのか？

私たちが地図上でデータを表示したい場合、投影を使用する必要があります。2Dの平面上に何かを表示するためには、3Dの地球上の座標を投影する必要があります。私たちがよく知っていて使用しているもっとも人気のある投影方法は、メルカトル投影（ウェブメルカトルまたはWGS84メルカトルとも呼ばれます）であり、最も人気のある座標系は世界測地系1984（WGS84）（またはEPSG:4326）です。この座標系は度に基づいており、経度は-180°から180°（西から東）、緯度は-90°から90°（南から北）までの範囲となっています。

• ポイントクラウド用のセグメント化ガイド「Segment Anything 3D for Point Clouds Complete Guide (SAM 3D)」
• 2024年に探すべき6つのリモートAIジョブ
• 「ウェアラブルデータによるコロナ感染予測」

2D平面への各投影にはいくつかの歪みがあります。メルカトル投影は角度が地球上のオブジェクト間で保持される保形マップ投影です。緯度が0°より上（または0°より下）に高くなるほど、領域と距離の歪みが大きくなります。ボロノイ図はシード間の距離に大きく依存しているため、図を生成する際に同じ歪みのエラーが転送されます。

地球は不規則な形を持つ楕円体ですが、私たちの目的のためには球体の形で近似できます。球体上でボロノイ図を生成することで、球体の表面上の弧に基づいて正確に距離を計算することができます。後で、生成された球状ポリゴンを投影された2D座標にマッピングすることができ、これらのセルを定義する2つのシードを結ぶ線に対して垂直な線がセルを分けることができることを確認することができます。

下記に、上記で説明した角度と距離の問題を示します。線は同じ点で交差しているにもかかわらず、ボロノイセルの形状と角度は異なります。

もう一つの問題は、2Dボロノイ図を使用した場合、世界の異なる地域（つまり、同じ緯度上にない地域）の領域を比較することができないことです。なぜなら、エリアが重度に歪むからです。

下記の例で使用されたコードと共にフルのJupyterノートブックはGitHubで見つけることができます。ここでは簡潔にするためにいくつかの関数が省略されています。

前提条件

必要なライブラリをインストールする

pip install -q srai[voronoi,osm,plotting] geodatasets

必要なモジュールと関数をインポートする

import geodatasetsimport geopandas as gpdimport matplotlib.pyplot as pltimport plotly.express as pxfrom shapely.geometry import MultiPoint, Pointfrom shapely.ops import voronoi_diagramfrom srai.regionalizers import VoronoiRegionalizer, geocode_to_region_gdf

最初の例

地球上には以下の6つのポイントを定義しましょう：北極と南極、そして赤道上に4つのポイントです。

earth_points_gdf = gpd.GeoDataFrame(    geometry=[        Point(0, 0),        Point(90, 0),        Point(180, 0),        Point(-90, 0),        Point(0, 90),        Point(0, -90),    ],    index=[1, 2, 3, 4, 5, 6],    crs="EPSG:4326",)

Shapelyライブラリを使用してvoronoi_diagramを使ってボロノイ図を生成する

def generate_flat_voronoi_diagram_regions(    seeds_gdf: gpd.GeoDataFrame,) -> gpd.GeoDataFrame:    points = MultiPoint(seeds_gdf.geometry.values)    # 2Dダイアグラムを生成    regions = voronoi_diagram(points)    # ジオメトリをGeoDataFrameにマップする    flat_voronoi_regions = gpd.GeoDataFrame(        geometry=list(regions.geoms),        crs="EPSG:4326",    )    # シードデータフレームからインデックスを適用する    flat_voronoi_regions.index = gpd.pd.Index(        flat_voronoi_regions.sjoin(seeds_gdf)["index_right"],        name="region_id",    )    # 地球の境界にクリップする    flat_voronoi_regions.geometry = flat_voronoi_regions.geometry.clip_by_rect(        xmin=-180, ymin=-90, xmax=180, ymax=90    )    return flat_voronoi_regions

earth_poles_flat_voronoi_regions = generate_flat_voronoi_diagram_regions(    earth_points_gdf)

sraiライブラリのVoronoiRegionalizerを使ってボロノイ図を生成する。内部では、scipyライブラリから実装されたSphericalVoronoiを使用し、WGS84座標を球面座標系に正しく変換します。

earth_points_spherical_voronoi_regions = VoronoiRegionalizer(    seeds=earth_points_gdf).transform()

両者のプロットでの違いを見てみましょう。

最初に目にすることは、2Dのヴォロノイ図は地球上を一周せず、フラットなデカルト座標平面上で機能するためです。球体のヴォロノイ図は地球全体を適切にカバーし、対子午線（経度が180°から-180°に切り替わる線）で切れなくなっています。

差を数値的に特定するために、ポリゴンの形状の差を測るためにIoU（Intersection over Union）メトリックス（またはJaccard指数）を計算することができます。このメトリックスの値は、0から1の間であり、0は重なりがなく、1は完全な重なりを意味します。

def calculate_iou(    flat_regions: gpd.GeoDataFrame, spherical_regions: gpd.GeoDataFrame) -> float:    total_intersections_area = 0    total_unions_area = 0    # すべての領域を繰り返す    for index in spherical_regions.index:        # 一致する球体とフラットなヴォロノイ領域を見つける        spherical_region_geometry = spherical_regions.loc[index].geometry        flat_region_geometry = flat_regions.loc[index].geometry        # 交差領域の面積を計算する        intersections_area = spherical_region_geometry.intersection(            flat_region_geometry        ).area        # 共用領域の面積を計算する        # 代替コード：        # spherical_region_geometry.union(flat_region_geometry).area        unions_area = (            spherical_region_geometry.area            + flat_region_geometry.area            - intersections_area        )        # 総合計に追加する        total_intersections_area += intersections_area        total_unions_area += unions_area    # 交差領域を合計領域で割る    return round(total_intersections_area / total_unions_area, 3)

calculate_iou(    earth_points_flat_voronoi_regions, earth_points_spherical_voronoi_regions)

計算された値は0.423であり、かなり低く、大規模なスケールでは、これらのポリゴンは互いに異なっていることが上のプロットで簡単に確認できます。

実データの例：AED（自動外部除細動器）位置を使用して地球を分割する

この例で使用されるデータはOpenAEDMapから取得され、OpenStreetMapのデータに基づいています。準備されたファイルには、重複したノードが重なって定義されていないポジション（正確には80694個）がフィルタリングされています。

# AEDの位置をGeoDataFrameにロードaed_world_gdf = gpd.read_file(    "https://raw.githubusercontent.com/RaczeQ/medium-articles/main/articles/spherical-geovoronoi/aed_world.geojson")

AEDのためのヴォロノイ図を生成します

aed_flat_voronoi_regions = generate_flat_voronoi_diagram_regions(aed_world_gdf)aed_spherical_voronoi_regions = VoronoiRegionalizer(    seeds=aed_world_gdf, max_meters_between_points=1_000).transform()

これらのヴォロノイ図を比較しましょう。

プロットを見ると、違いはかなり明らかです。2Dバージョンではすべての境界線が直線ですが、球面のものはWGS84座標では曲がって見えます。また、フラットなバージョンでは、たくさんの領域が極地に収束しているのが明確にわかります（正射投影は南極に焦点を当てています）、一方、球面のものでは収束しません。また、最初の例で言及された対子午線周りの連続性の違いも見られます。ニュージーランドから出ている領域はフラットなバージョンでは突然切れています。

IoU値を見てみましょう：

calculate_iou(aed_flat_voronoi_regions, aed_spherical_voronoi_regions)

計算された値は0.511で、最初の例よりわずかに良くなっていますが、ポリゴンの一致率はおおよそ50％です。

都市スケールでの拡大

より小さなスケールでの違いを見てみましょう。ロンドンにあるすべてのAEDを選択し、プロットします。

greater_london_area = geocode_to_region_gdf("Greater London")aeds_in_london = aed_world_gdf.sjoin(greater_london_area)

calculate_iou(    aed_flat_voronoi_regions.loc[aeds_in_london.index],    aed_spherical_voronoi_regions.loc[aeds_in_london.index],)

値は0.675です。改善されていますが、まだ顕著な違いがあります。AEDがより密に配置されているため、形状と距離が縮小され、射影された2D平面と球面で計算されたボロノイ図の違いが小さくなります。

いくつかの個別の例を重ね合わせて見てみましょう。

これらのポリゴンの領域はほとんど一致していますが、角度や形状の違いが見られます。これらの相違点は空間分析において重要であり、結果を変える可能性があります。対象領域が大きくなるほど、違いも大きくなります。

まとめ

平面のボロノイ図よりも、球面のボロノイ図が地理空間ドメインでの使用に適している理由がわかったと思います。

現在、このドメインでのほとんどの分析は、射影された2D平面でのボロノイ図を使用して行われており、誤った結果につながる可能性があります。

長い間、Pythonで地理空間データサイエンティストやアナリストが使用できる簡単な球面のボロノイ図の解決策はありませんでした。今は、1つのライブラリをインストールするだけで簡単になりました。2Dの解決策よりも計算時間は少し長くなりますが、球面座標への射影と射影中のポリゴンの適切なクリッピングが必要です。ただし、分析結果の精度を保持したい場合、これは重要ではありません。JavaScriptユーザーには、すでに球面のボロノイ図のD3.js実装が利用可能です。

免責事項

私はsraiライブラリのメンテナの一人です。

We will continue to update VoAGI; if you have any questions or suggestions, please contact us!