非GISデータベースでGIS計算を実行する方法

導入：

場合によっては、ユーザーの位置のマッピングや地理データの分析など、アプリケーション内で地理空間機能を実行する必要に迫られることがあります。これらのタスクには、GDAL、Shapely、Python 用の Geopandas など、言語固有のライブラリが多数用意されています。

あるいは、地理空間機能はデータベースを通じて実装することもできます。たとえば、PostgreSQL などのリレーショナル データベースで PostGIS 拡張機能を使用したり、Azure CosmosDB などの分散データベースで空間データ型のネイティブ サポートを活用したりすることができます。

ただし、Redis や Google Spanner などのバックエンド ストレージが空間クエリをネイティブにサポートしておらず、大規模な地理空間クエリを処理する必要がある場合は、この記事が役立ちます。

選択肢は何ですか?

空間データを処理する別のマイクロサービスをいつでも構築できますが、このオプションでは、追加のアプリケーションを維持するオーバーヘッドが発生することがよくあります。別のアプローチは、S2 や H3 などの地理インデックス ライブラリを使用することです。Google が開発した S2 はヒルベルト曲線に基づいており、Uber が開発した H3 は測地線離散グローバル グリッド システムに基づいています。S2 と H3 には多くの類似点があります。どちらも特定の領域をセルに分割し、64 ビットの整数を使用してこれらのセルをインデックスします。

ただし、主な違いはセルの形状にあります。S2 は正方形のセルを使用するのに対し、H3 は六角形のセルを使用します。一部のアプリケーションでは、H3 の方がパフォーマンスが優れている場合があります。ただし、全体的にはどちらのライブラリでも十分です。この記事では S2 を使用しますが、H3 を使用しても同様の機能を実行できます。

Google S2 ライブラリの基本概念

• セル: S2 は球体をセルに分割し、各セルには一意の 64 ビット識別子が割り当てられます。

• セル レベル: 階層により、大きな領域から小さな正確な領域まで、さまざまなレベルの詳細が可能になります。各レベルは異なる解像度を表します。

  
  

  • レベル 0: 地球の表面の大部分を覆う最大のセル。

    
    

  • 上位レベル: セルは段階的に小さな象限に分割されます。たとえば、レベル 1 のセルはそれぞれ 4 つのレベル 2 のセルに分割されます。

    
    

  • 解像度と領域: レベルが高いほど、解像度が高くなり、セル領域が小さくなります。この階層により、さまざまな詳細レベルで正確なインデックス作成とクエリが可能になります。

    
    

以下の表は、さまざまなセル レベルとそれに対応する領域を示しています。

提供された表から、S2 を使用すると 0.44 cm^2 までのマッピング精度を実現できることは明らかです。S2 セルの各四角形内には、同じ親を共有する子セルが存在し、階層構造を示しています。セルのレベルは静的な値 (すべてのセルに同じレベルを適用) にすることも、S2 が最適な解像度を決定する動的な値にすることもできます。

最近傍点の計算

まずは例から始めましょう。シアトル地域に近接サービスのような機能を提供するアプリケーションを作成するとします。指定された範囲内のコーヒーショップのリストを返します。これらの操作を実行するには、このタスクを 4 つのサブタスクに分割します。

• シアトルの地図を読み込んでいます
• シアトルの地図上で S2 セルを視覚化する
• いくつかのコーヒーショップの場所をデータベースに保存する
• 最寄りのコーヒーショップを検索

シアトルマップを読み込んでいます

Google マップを読み込むには、gmplot ライブラリを使用します。このライブラリを読み込むには、Google マップ API キーが必要です。API キーを生成するには、ここ の手順に従ってください。

 import gmplot import const # plot seattle with zoom level 13 gmap = gmplot.GoogleMapPlotter(47.6061, -122.3328, 13, apikey=const.API_KEY) # Draw the map to an HTML file: gmap.draw('map.html')

上記のコードは、以下に示すような map.html ファイルを生成します。

シアトルマップ上で S2 セルを視覚化する

マップができたので、マップ用の S2 セルをいくつか描画してみましょう。

 from s2 import * import gmplot # plot seattle with zoom level 13 gmap = gmplot.GoogleMapPlotter(47.6061, -122.3328, 13, apikey=const.API_KEY) areatobeplotted = [ (47.64395531736767,-122.43597221319135), (47.51369277846956,-122.43597221319135), (47.51369277846956,-122.24156866779164), (47.64395531736767,-122.24156866779164), (47.64395531736767,-122.43597221319135) ] region_rect = S2LatLngRect( S2LatLng.FromDegrees(47.51369277846956,-122.43597221319135), S2LatLng.FromDegrees(47.64395531736767, -122.24156866779164)) coverer = S2RegionCoverer() coverer.set_min_level(8) coverer.set_max_level(15) covering = coverer.GetCovering(region_rect) geoms = 0 for cellid in covering: new_cell = S2Cell(cellid) vertices = [] for i in range(0, 4): vertex = new_cell.GetVertex(i) latlng = S2LatLng(vertex) vertices.append((latlng.lat().degrees(), latlng.lng().degrees())) gmap.polygon(*zip(*vertices), face_color='pink', edge_color='cornflowerblue', edge_width=5) geoms+=1 gmap.polygon(*zip(*areatobeplotted), face_color='red', edge_color='green', edge_width=5) print(f"Total Geometries: {geoms}") gmap.draw('/tmp/map.html')

 Output: Total Geometries: 273

上記のコードでは、まず Google マップ プロッターをシアトル周辺を中心に配置します。 S2RegionCovererでは、リージョン カバラーを初期化して、最小レベル 8 から最大レベル 15 までの動的レベルを設定します。これにより、S2 はすべてのセルを特定のセル サイズに動的に適合させ、最適なフィットを実現できます。 GetCoveringメソッドは、シアトル周辺の四角形のカバリングを返します。

次に、各セルを反復処理して、セルの頂点を計算し、マップ上にプロットします。生成されるセルの数は約 273 個に抑えます。最後に、入力された四角形を赤でプロットします。このコードは、以下に示すように、シアトル マップの/tmp/map.htmlに S2 セルをプロットします。

いくつかのコーヒーショップの場所をデータベースに保存する

コーヒーショップのデータベースを S2 セル識別子とともに生成してみましょう。これらのセルは任意のデータベースに保存できます。このチュートリアルでは、SQLite データ データベースを使用します。以下のコード サンプルでは、SQLite データベースに接続して、 Id 、 name 、 cell_idの 3 つのフィールドを持つCoffeeShopsテーブルを作成します。

前の例と同様に、 S2RegionCovererを使用してセルを計算しますが、今回はポイントをプロットするための固定レベルを使用します。最後に、計算された ID が文字列に変換され、データベースに保存されます。

 import sqlite3 from s2 import S2CellId,S2LatLng,S2RegionCoverer # Connect to SQLite database conn = sqlite3.connect('/tmp/sqlite_cells.db') cursor = conn.cursor() # Create a table to store cell IDs cursor.execute('''CREATE TABLE IF NOT EXISTS CoffeeShops ( id INTEGER PRIMARY KEY AUTOINCREMENT, name TEXT, cell_id TEXT )''') coverer = S2RegionCoverer() # Function to generate S2 cell ID for a given latitude and longitude def generate_cell_id(latitude, longitude, level=16): cell=S2CellId(S2LatLng.FromDegrees(latitude, longitude)) return str(cell.parent(level)) # Function to insert cell IDs into the database def insert_cell_ids(name,lat,lng): cell_id = generate_cell_id(lat, lng) cursor.execute("INSERT INTO CoffeeShops (name, cell_id) VALUES (?, ?)", (name, cell_id)) conn.commit() # Insert cell IDs into the database insert_cell_ids("Overcast Coffee", 47.616656277302155, -122.31156460382837) insert_cell_ids("Seattle Sunshine", 47.67366852914391, -122.29051997415843) insert_cell_ids("Sip House", 47.6682364706238, -122.31328618043693) insert_cell_ids("Victoria Coffee",47.624408595334536, -122.3117362652041) # Close connection conn.close()

この時点で、セル レベルで選択された解像度によって決定されるセル ID とともにコーヒー ショップを保存するデータベースができました。

最寄りのコーヒーショップを検索

最後に、大学地区エリアのコーヒーショップを検索してみましょう。

 import sqlite3 from s2 import S2RegionCoverer,S2LatLngRect, S2LatLng # Connect to SQLite database conn = sqlite3.connect('/tmp/sqlite_cells.db') cursor = conn.cursor() # Function to query database for cells intersecting with the given polygon def query_intersecting_cells(start_x,start_y,end_x,end_y): # Create S2RegionCoverer region_rect = S2LatLngRect( S2LatLng.FromDegrees(start_x,start_y), S2LatLng.FromDegrees(end_x,end_y)) coverer = S2RegionCoverer() coverer.set_min_level(8) coverer.set_max_level(15) covering = coverer.GetCovering(region_rect) # Query for intersecting cells intersecting_cells = set() for cell_id in covering: cursor.execute("SELECT name FROM CoffeeShops WHERE cell_id >= ? and cell_id<=?", (str(cell_id.range_min()),str(cell_id.range_max()),)) intersecting_cells.update(cursor.fetchall()) return intersecting_cells # Query for intersecting cells intersecting_cells = query_intersecting_cells(47.6527847,-122.3286438,47.6782181, -122.2797203) # Print intersecting cells print("Intersecting cells:") for cell_id in intersecting_cells: print(cell_id[0]) # Close connection conn.close()

 Output: Intersecting cells: Sip House Seattle Sunshine

以下はセルの視覚的表現です。簡潔さを保つために、以下の視覚化コードは追加されていません。

すべての子セルと親セルはプレフィックスを共有しているので、最小値と最大値の間のセル範囲をクエリして、その2つの値の間にあるすべてのセルを取得できます。例では、同じ原理を使用してコーヒーショップをクエリします。

結論：

この記事では、地理インデックスを使用して、地理空間クエリをサポートしていないデータベースに地理空間データを保存およびクエリする方法を説明しました。これは、2 点間のルートを計算したり、最も近い近傍を取得したりするなど、複数のユース ケースにさらに拡張できます。

通常、地理インデックス付きデータベースのクエリでは、データに対して追加の後処理を実行する必要があります。ノードに過負荷をかけないように、クエリと後処理のロジックを慎重に検討する必要があります。

参考文献:

• Google マップ プロッター - https://github.com/gmplot/gmplot/wiki/GoogleMapPlotter
• S2 開発者ガイド - http://s2geometry.io/devguide/
• SQLite - https://docs.python.org/3/library/sqlite3.html
• Azure Cosmos DB 地理空間 - https://learn.microsoft.com/en-us/azure/cosmos-db/nosql/query/geospatial?tabs=javascript
• GDAL - https://gdal.org/index.html
• Shapley - https://shapely.readthedocs.io/en/stable/manual.html
• ジオパンダ - https://geopandas.org
• ポスギス - https://postgis.net/
• Google スパナー - https://cloud.google.com/spanner?hl=en
• レディス - https://redis.io/