PostGIS拡張を導入したAlloy DBにGCEからAuth Proxyを介して接続する

はじめに

こんにちは、GIクラウドのスペランカーYoshioです。
今回はAlloy DBと銘打ちつつ、PostGIS拡張で遊ぶ準備をしているのが実態の記事です。さすがにAlloy DBっぽさ皆無なのはよろしくないので、Alloy DB Auth Proxyを使ってお茶を濁しています。

点を打つくらいなら簡単な地理情報ですが、多角形に対する内外判定とかやりだすと途端に難易度が上がります。地球は回転楕円体なので、そこそこの精度で経緯度から2点間の距離を算出することですら意外と難しいものです。算数感覚でポリゴンの面積が欲しいと言われて「地球は楕円体なんだぞ！」って心の中で絶叫した人は僕だけじゃないはず。まぁ、実装する楽しさはあるのですが、そんなことに工数を割けるPJなんてそうはありません。

そこで、PostgreSQLのPostGIS拡張の出番です。この拡張機能を使って地理情報の計算処理を簡単にできるようにしましょうというのがこの記事の趣旨です。

ゴール

• GCEからAuth Proxyを介したプライマリインスタンスへの接続で LINESTRING 形式のデータを投入する
• GCEからAuth Proxyを介したリードプールインスタンスへの接続で LINESTRING の長さを取得

※ LINESTRING はPostGISでの折線のことです。点と点で接続されています。

サマリ

以下、この記事の目次です。僕が参考にした公式ドキュメントも併記しておきます。

1. Alloy DB Auth Proxy導入
   • About the AlloyDB Auth proxy
   • Connect using the AlloyDB Auth proxy
   • Best practices for using the AlloyDB Auth proxy
2. PostGIS拡張の導入とテーブルの作成
   • Supported database extensions
3. LINESTRINGデータの投入
4. データの読み込みと距離の算出

前提条件

Alloy DBの環境構築については。弊社ブログAlloy DB試してみたをご参照ください。

---

Alloy DB Auth Proxy導入

About the AlloyDB Auth proxyによると、Alloy DB Auth Proxyは”Alloy DBへの安全な接続を提供する機能”です。

以下のような前提条件によると、同一VPCネットワーク内で完結するこのサンプルではAuth Proxyはそんなに必要性がなさそうなのですが、条件もそろっていますし、せっかく提供されているので使ってみます。

It must have network visibility to the VPC network where the instances you want to connect to reside. （接続先のインスタンスが存在するVPCネットワークへのネットワーク可視性が必要です。）

Best practices for using the AlloyDB Auth proxyを参考にしつつ、公式ドキュメントのConnect using the AlloyDB Auth proxyに従って進めます。手順は以下のようになります。

1. GCEデフォルトのサービスアカウントに接続用の鍵ファイルを追加
2. Auth ProxyクライアントのDL
3. 接続

GCEデフォルトのサービスアカウントにAlloy DB接続用の鍵ファイルを追加

ナビゲーションメニューの IAMと管理 から サービスアカウント ページを表示。 {SERVICE_PROJECT_NUMBER}-compute@developer.gserviceaccount.com のアクションメニューから 鍵を管理 を選択して鍵管理ページを表示。鍵を追加 プルダウンボタンがあるので、そこから 新しい鍵を作成 を選択します。追加する鍵はJSON形式です。
※ 取得した鍵ファイルは大切に保管してください。

任意の方法で鍵ファイルをAlloy DBに接続するGCEインスタンスにアップロードしてください。

Auth ProxyクライアントのDL

SSHでAlloy DBに接続するGCEインスタンスに接続してAuth Proxyをダウンロードし、実行権限を付与します。ちなみに検証で利用しているGCEインスタンスのOSは Debian GNU/Linux 11 です。

$ curl -o alloydb-auth-proxy https://storage.googleapis.com/alloydb-auth-proxy/v0.2.0/alloydb-auth-proxy.linux.amd64

$ chmod +x alloydb-auth-proxy

接続

Cloud Shell で対象となるAlloy DBの接続URIを取得します。予め設定してあるプライマリインスタンスとリードインスタンスのURIを取得します。
※ projects/{PROJECT_ID}/locations/{REGION_ID}/clusters/{CLUSTER_ID}/instances/{INSTANCE_ID} で作成可能です。

$ gcloud beta alloydb instances list --cluster={CLUSTER_ID} --region={REGION_ID} --project={PROJECT_ID}

NAME: projects/{PROJECT_ID}/locations/{REGION_ID}/clusters/{CLUSTER_ID}/instances/{READ_POOL_INSTANCE_ID}
INSTANCE_TYPE: READ_POOL
STATUS: READY

NAME: projects/{PROJECT_ID}/locations/{REGION_ID}/clusters/{CLUSTER_ID}/instances/{PRIMARY_INSTANCE_ID}
INSTANCE_TYPE: PRIMARY
STATUS: READY

Auth Proxyクライアントを起動します。今回は取得したプライマリインスタンスの接続Portを5432に、リードプールインスタンスの接続Portを5433に設定しています。

この例では使っていませんが、&をつけてバックグラウンド起動をするケースも多そうです。

$ ./alloydb-auth-proxy "projects/{PROJECT_ID}/locations/{REGION_ID}/clusters/{CLUSTER_ID}/instances/{PRIMARY_INSTANCE_ID}:5432" "projects/{PROJECT_ID}/locations/{REGION_ID}/clusters/{CLUSTER_ID}/instances/{READ_POOL_INSTANCE_ID}:5433" --credentials-file={鍵ファイルのパス}

起動に成功したら、以下のようなメッセージが表示されます。

[1] 3765 # PID
$ Authorizing with the credentials file at "{xxxx.json}"
[projects/{PROJECT_ID}/locations/{REGION_ID}/clusters/{CLUSTER_ID}/instances/{PRIMARY_INSTANCE_ID}:5432] Listening on 127.0.0.1:5432
[projects/{PROJECT_ID}/locations/{REGION_ID}/clusters/{CLUSTER_ID}/instances/{READ_POOL_INSTANCE_ID}:5433] Listening on 127.0.0.1:5433
The proxy has started successfully and is ready for new connections!

実際にプライマリインスタンスとリードプールインスタンスのそれぞれへの接続を確認します。以下はリードプールインスタンスへの接続例です。

$ psql -h 127.0.0.1 -p 5433 -U postgres

PostGIS拡張の導入とテーブルの作成

プライマリインスタンスにアクセスします。この場合、Portが5432となっています。

$ psql -h 127.0.0.1 -p 5432 -U postgres

1. sampleデータベースを切る
2. 作成したsampleデータベースに切り替え
3. sampleデータベースにPostGIS拡張を導入
4. geography型のカラムを持ったテーブルを作成
   • この例ではWGS84のLINESTRINGを持たせています

postgres=> CREATE DATABASE sample;
CREATE DATABASE
postgres=> \c sample;
sample=> CREATE EXTENSION IF NOT EXISTS postgis;
CREATE EXTENSION
sample=> CREATE TABLE routes(route_id serial PRIMARY KEY, route_path geography(LINESTRING, 4326));
CREATE TABLE

以下のようなエラーが表示されることがありますが、これは \c コマンドでテーブル切り替えをする前にPostGIS拡張を適用した場合に発生するようです。

ERROR:  type "geography" does not exist
LINE 1:         route_path GEOGRAPHY(LINESTRING, 4326),

---

LINESTRINGデータの投入

ちょうど手元にあったGarminで計測したトラッキングデータを使ってデータをLINESTRINGデータを作成します。GPXファイル形式でトラッキングデータをエクスポートし、そこから経緯度情報を取得します。
run.gpx というGPXファイルを読み込んでINSERTする例です。Prepared Statement使っていないのはご愛敬ということで。

Auth Proxyで接続しているので、接続先が 127.0.0.1 になっています。このプログラムをGCEで実行します。

import os
import math
import xml.etree.ElementTree as ET
from functools import reduce
import sqlalchemy

db = None

def store_route(route:list):
    query = convert_route_to_linestringquery(route)
    with db.connect() as conn:
        record = conn.execute(query)

def connect_tcp_socket() -> sqlalchemy.engine.base.Engine:
    pool = sqlalchemy.create_engine(
        sqlalchemy.engine.url.URL.create(
            drivername="postgresql+pg8000",
            username="postgres",
            password="**********",
            host="127.0.0.1",
            port=5432,
            database="sample",
        ),
        pool_size=5,
        max_overflow=2,
        pool_timeout=5,
        pool_recycle=60,
    )
    return pool

def init_db() -> sqlalchemy.engine.base.Engine:
    global db
    db = connect_tcp_socket()

def find_all_treking_latlon(path:str) -> list:
    with open(path, encoding="utf8") as gpx:
        tree = ET.parse(gpx)
        root = tree.getroot()
        return [(elem.attrib["lat"], elem.attrib["lon"]) for elem in root.iter() if elem.tag=="{http://www.topografix.com/GPX/1/1}trkpt"]

def convert_route_to_linestringquery(route:list) -> str:
    points = ",".join([pt[1]+" "+pt[0] for pt in route])
    return "INSERT INTO routes(route_path) VALUES (ST_GeographyFromText('SRID=4326;LINESTRING("+points+")'))"

# 検算用の経緯度から距離を算出する検算用プログラム
POLE_RADIUS = 6356752    # 極半径(短半径)
EQUATOR_RADIUS = 6378137 # 赤道半径(長半径)
E = 0.081819191042815790 # 離心率
E2= 0.006694380022900788 # 離心率の２乗

class Coordinate:
    def __init__(self, latitude, longitude):
        self.latitude  = float(latitude)
        self.longitude = float(longitude)

def calc_distance(point_a, point_b):
    a_rad_lat = math.radians(point_a.latitude)
    a_rad_lon = math.radians(point_a.longitude)
    b_rad_lat = math.radians(point_b.latitude)
    b_rad_lon = math.radians(point_b.longitude)
    m_lat = (a_rad_lat + b_rad_lat) / 2 # 平均緯度
    d_lat = a_rad_lat - b_rad_lat # 緯度差
    d_lon = a_rad_lon - b_rad_lon # 経度差
    W = math.sqrt(1-E2*math.pow(math.sin(m_lat),2))
    M = EQUATOR_RADIUS*(1-E2) / math.pow(W, 3) # 子午線曲率半径
    N = EQUATOR_RADIUS / W # 卯酉線曲率半径
    d = math.sqrt(math.pow(M*d_lat,2) + math.pow(N*d_lon*math.cos(m_lat),2))
    return d

def calc_route_distance(route:list) -> float:
    coordinates = [Coordinate(pt[0], pt[1]) for pt in route]
    l = len(coordinates)
    distance = 0
    for i in range(l):
        if i+1 == l:
            break
        distance += calc_distance(coordinates[i], coordinates[i+1])
    return distance

if __name__ == "__main__":
    init_db()
    route = find_all_treking_latlon("run.gpx")
    store_route(route)
    print("route distance: " + str(calc_route_distance(route)))

ついでに、トラッキングデータの距離の検算で用いる値の算出処理も実行しています。ここではヒュベニの公式を使っています。ソースコードはヒュベニの公式（Hybeny’s Distance Formula）による緯度・経度の2点間の距離の算出からお借りしています。
この例では、9012.172484353374 メートルという距離が算出されました。参考元のサイトによるとヒュベニの公式自体は近似値を算出するものらしいので多少の誤差が出るでしょう。

---

データの読み込みと距離の算出

登録したデータをリードプールから読み込んでみます。GCEから psql コマンドでリードプールにアクセスします。

$ psql -h 127.0.0.1 -p 5433 -U postgres

ST_Length 関数を使ってLINESTRINGの長さを求めてみます。

sample=> select st_length(route_path,false) from routes;
     st_length
-------------------
 9024.240971095775
(1 row)

sample=> select st_length(route_path,true) from routes;
     st_length
-------------------
 9012.172484467554
(1 row)

ST_Length 関数の第2引数のフラグは以下のような意味を持ちます。

フラグ	説明
false	球面計算のため精度が低い
true（or 未指定）	回転楕円体面モデルでの計算のため精度が高い

先ほど取得した検算用の値 9012.172484353374 と比較してもおおよそ問題ない数値を得ることができました。

---

おわりに

都度データ登録する必要がない場合も、以下のようなPostGISのST_Distanceサンプルのように演算結果だけ取得する方法も可能です。PostGISと読み取りプール インスタンスとの組み合わせ方次第では、地理情報処理の専門知識がなくてもアプリケーション作成の幅を広げることができます。急にポリゴンの面積を出して欲しいと言われてももう大丈夫。

sample=> SELECT ST_Distance(gg1, gg2) As spheroid_dist, ST_Distance(gg1, gg2, false) As sphere_dist FROM (SELECT 'SRID=4326;POINT(-72.1235 42.3521)'::geography as gg1, 'SRID=4326;LINESTRING(-72.1260 42.45, -72.123 42.1546)'::geography as gg2) As foo;
 spheroid_dist | sphere_dist
---------------+--------------
  123.80207675 | 123.47573692
(1 row)

前回記事のCloud Runからの呼び出しを活用すれば、さらに地理情報をアプリケーションとしての使い道も広がってくるのではないでしょうか。まぁ、Alloy DBでやる必要はあるかと問われると困ってしまうのですが。

このような感じで、途中から完全にPostGISの話になってしまいましたが、全部DBで管理したい。スケーラビリティも欲しい。地理情報の処理もしたいという欲張りなお客さんにうってつけなAlloy DB、おススメです。

---

GI Cloud は事業の拡大に向けて一緒に夢を追う仲間を募集しています

当社は「クラウドで日本のIT業界を変革し、世の中をもっとハッピーに」をミッションに掲げ、Google Cloudに特化した技術者集団として、お客様にコンサルティングからシステム開発、運用・保守まで一気通貫でサービスを提供しています。

まだ小規模な事業体ですが、スタートアップならではの活気と成長性に加えて、大手総合商社である伊藤忠グループの一員としてやりがいのある案件にもどんどんチャレンジできる環境が整っています。成長意欲の高い仲間と共にスキルを磨きながら、クラウドの力で世の中をもっとハッピーにしたい。そんな我々の想いに共感できる方のエントリーをお待ちしています。

採用ページ