RDBMSで位置情報を扱おう ～PostGIS超入門 2024/10/07 有限会社アートライ／日本仮想化技術株式会社 坂井 恵 第4回 爆速DB「PG-Strom」勉強会 https://bakusokudb.connpass.com/event/330639/

自己紹介 • 坂井 恵 @sakaik • 有限会社アートライ • 日本仮想化技術株式会社「爆速DB Powered by PG-Strom」チームに参画中 ■GIS関連の略歴 • MySQL 8.0(2018)で地理情報データを扱う機能が大幅に充実されたことをきっ かけに、地理情報データに興味 • 「足で稼ぐ地理情報」がモットー • 測量士補

今日のテーマとターゲット ■タイトル RDBMSで位置情報を扱おう～PostGIS超入門 ■ターゲット • RDBMSで位置情報（地理情報）を扱うことに興味がある人 • RDBMSの操作はある程度経験があることを前提 • これまでRDBMSで地理情報を扱ったことがない人 をメインターゲットとして、入門エッセンスをお届けします ※今日はPG-Stromの話はほぼありません

今日の構成 • 地理情報入門 • RDBMS(PostGIS)で地理情報を扱う方法の例 ぎゅっと概要をお伝えするために、多少厳密さを欠く部分があります。 ご了承ください。本質には影響がないようお伝えするつもりです。

地理情報入門～地理データとは？ • 地理情報データが表すもの • 地球上の場所（POINT) • 全国の某ファストフード店舗の場所情報 • 公衆トイレ情報 • 人がどこにいるか（人流データなど） • 地球上の点をいくつか繋いだ線(LINESTRING) • 河川 • 鉄道 • 地球上にあるエリア(POLYGON) • 自治体（市区町村） • 湖沼 • 建物の形や敷地の形 ※必ずしも「地球上」である必要はない のですが、今日は話をシンプルにするた めによくつかう対象として「地球上」と 表現しました。

「POINT」がすべての基本 • POINT • (緯度＋経度)などで表す • 詳細はこのあと説明します • （他の方法もあります） • このページではいくつかの点を仮に、A,B,C,D,Eと呼びます • LINESTRING • 複数のPOINTを繋いだもの。 • 例）A-B-C-D, B-D など • POLYGON • 複数のPOINTを繋いでエリアを示したもの • 例) A-B-C-D-A, C-D-E-C など • 中をくり抜く仕組みもある

「POINT」がすべての基本 • POINTを知るための基礎知識 • 地球の形と大きさ • 様々な「地球のモデル」 • 地理座標と投影座標 ※今日は詳細の説明はせず「こういうのがあります」と紹介するに留めま すが、定義に関する重要な話なので、今日出てきた話は何かのタイミング でぜひ調べてみてください。

地球の形と大きさ • 地球はデコボコ • 「回転楕円体」としてモデル化する • 長半径、扁平率 の2つの数値により定義できる • モデル＝測地系 • 日本でよく使われる測地系 • JGD2011 (SRID=6668) • WGS84 (SRID=4326) • 同じ緯度経度数値でも、測地系によって指し示す場所が異なる

たとえばなし：「温度型」 • DBの型として仮に「温度型」というのがあると想像してください。 • 単なる数値ですが、便利な関数があるとかそんなんで独自の「型」が用意さ れているとか • いま、温度型のカラムに「36.74」と入れたとします • さて、この列の値は何度を示していますか？ • → 摂氏なのか華氏なのかの情報がこの型には必要ですね！ • 例）INSERT INTO ..... VALUES ( MkTmprtr(36.74, 'C') ); • 例）SELECT Tmpr_AsUnit( col1, 'F' ) FROM ...; → 98.132 • 地理情報も同じです！ • 緯度経度だけでなく、それが何の測地系なのかの情報が常に必要です

緯度と経度 • 日本で地理院などが配布している地理データの測地系 • 大概は JGD2011 • 長半径 6,378,137m、扁平率 1/298.257222101 • GPSなど世界的にはWGS84が使われることが多い • 長半径 6,378,137m、扁平率 1/298.257223563 →実務上は（センチ／ミリ単位の厳密な議論を行うのでなければ）ほぼ 一緒と考えて良い （個人的には「だいたい一緒なのでテキトーでいいや」ではなく「だい たい一緒だからこそ正しく使おう」とお伝えしたいです） ※データとして扱うときには、測地系の名前ではなく、測地系ごとに割り振られたSRID(SRS_ID) の数字を使います

緯度経度の「度」いろいろ • データとして扱うときは「度」で扱うことが多い • 度分の例 → 北緯 35度 39分 20.5482秒 東経 139度 42分 0.7158秒 北緯 35.6557082 度 東経139.7001994 度 北緯 35度 39.34247分 東経 139度 42.01193分 ＝ ＝ 度分秒 度分 度

(話題) 地理座標と投影座標 今日は詳細はお話しませんが、 位置を緯度経度で表す方法の他に、メートルで表す方法もあります（投影 座標）。 日本では「(日本の)平面直角座標系」がよく使われます（オープンなデータ にも平面直角座標系の値で公開されているものがあります） ※個人的に、この平面直角座標系の仕組みが大好きです。ぜひ機会あれば触れてみてください

ここまでのまとめ • 地球上の位置を表すために「緯度、経度」を使う • 緯度経度には「どの測地系で」という情報が必要 • 各測地系には SRID(SRSID)というIDが振られている • 日本でよく使うのは JGD2011(6668)、WGS84(4326) • 緯度経度はシステムでは「度」の単位で扱うことが多い

PostgeSQLでGIS機能を使う(準備) PostGISというextensionを入れる必要があります。 aptやyumで PostgreSQLとPostGISをインストール後、 PostgreSQLのスキーマ内で postgres=# create extension postgis;

テーブル作成とデータ登録 • GEOMETRY型の列を作る • ジオメトリ型に変換しながら登録 PostGISでJGD2011の時は 経度、緯度の順 経度緯度の間にコンマは使わない(ホワイトスペース) 測地系（SRID）を指定する postgres=# CREATE TABLE t1 (id integer, g geometry); postgres=# INSERT INTO t1 VALUES (1, ST_GeomFromText('POINT(135 35)', 6668)); INSERT 0 1 GEOMETRY型のほかに、 より限定した型やSRIDを指 定する方法もあります。 今日はこれで。

緯度経度の順番に注意 • PostGIS(JGD2011)では経度、緯度の順 • 逆に指定してもエラーにならないので注意 • MySQLの場合は同様の前提で「緯度、経度」の順 • MySQLの場合は緯度経度を逆に指定するとエラーにしてくれる mysql> INSERT INTO g1 VALUES (1, ST_GeomFromText('POINT(135 35)', 6668)); ERROR 3617 (22S03): Latitude 135.000000 is out of range in function st_geomfromtext. It must be within [-90.000000, 90.000000]. INSERT INTO t1 VALUES (2, ST_GeomFromText('POINT(35 135)', 6668)); INSERT INTO t1 VALUES (2, ST_GeomFromText('POINT(35 135)', 6668)); 緯度経度が逆！ ※正確に言うと逆にしてエラーになるわけではなく範囲を超えたらエラー。日本近郊に於いては逆にすると必ず緯度が範囲を超える。

地理データの参照 • 内部ではバイナリ形式で格納されている • 人間の目で理解しやすい形にする「ST_AsText()」覚えておくと便利 • 併せて、その値のSRIDを知る「ST_SRID()」も test=# SELECT id, g, ST_AsText(g), ST_SRID(g) FROM t1; id | g | st_astext | st_srid ----+----------------------------------------------------+---------------+--------- 1 | 01010000200C1A00000000000000E060400000000000804140 | POINT(135 35) | 6668 2 | 01010000200C1A000000000000008041400000000000E06040 | POINT(35 135) | 6668

参考：POINTデータもうひとつの登録方法 • 緯度と経度の値を持っている場合は、以下のように MakePoint してからSRIDを与える方法もあります。 INSERT INTO t1 VALUES (5, ST_SetSRID(ST_MakePoint(134.95, 35.02), 6668));

LINEとPOL YGONの登録方法 • LineStringの登録方法 • Polygonの登録方法 カッコが二重になっていることに注目 INSERT INTO t1 VALUES (4, ST_GeomFromText('Polygon((135.0 35.0, 135.01 35.0, 135.01 35.007, 135.005 35.015, 135.0 35.01, 135.0 35.0))',6668)); INSERT INTO t1 VALUES (3, ST_GeomFromText('LineString( 140.09475 35.8847, 140.12239 35.8719, 140.14917 35.84394, 140.21508 35.8623, 140.23534 35.85507)',6668));

ここまでのデータの確認 • ST_AsText()を使いこなしましょう test=# SELECT id, ST_AsText(g), ST_SRID(g) FROM t1 ORDER BY id; id | st_astext | st_srid ----+---------------------------------------------------------------------------------------------------------+--------- 1 | POINT(135 35) | 6668 2 | POINT(35 135) | 6668 3 | LINESTRING(140.09475 35.8847,140.12239 35.8719,140.14917 35.84394,140.21508 35.8623,140.23534 35.85507) | 6668 4 | POLYGON((135 35,135.01 35,135.01 35.007,135.005 35.015,135 35.01,135 35)) | 6668 (4 rows)

よく使うかもしれない情報源(1) • 国土数値情報ダウンロードサイト • https://nlftp.mlit.go.jp/ksj/ • 自治体ポリゴンや地域・地形等のデータ が多く公開されている

よく使うかもしれない情報源(2) • e-Stat 政府統計の総合窓口 • https://www.e-stat.go.jp/ • 人口や経済センサス、その他各種統計 調査の結果などが、主に自治体単位の 数値として公開されている • メッシュPolygonデータ(3次～5次)も ここに

地理情報データの活用方法の例 • 指定地点から周囲3500m以内にあるPOINTを抽出 距離を「度」で与える必要があるのが直感に反する 111319.5は地球の半径から計算した1度あたりの距離 SELECT id, ST_AsText(g) FROM t1 WHERE ST_DWithin(g, ST_SetSRID(ST_MakePoint(134.97,34.999),6668), 3500/111319.5); id | st_astext ----+--------------------------------------------------------------------------- 1 | POINT(135 35) 4 | POLYGON((135 35,135.01 35,135.01 35.007,135.005 35.015,135 35.01,135 35)) 5 | POINT(134.95 35.02) test=# SELECT 6378137*2*3.14159265358979/360; ?column? ----------------------- 111319.49079327345790

地理情報データの活用方法の例 • 指定したPolygon内にあるPOINTを抽出 SELECT p.id, ST_AsText(p.g) FROM points p JOIN area a ON (ST_Within(p.g, a.g)) WHERE a.name='渋谷区';

地理情報データの活用方法の例 • ST_Distance()を使ってジオメトリ間の距離を算出する方法も ある GEOGRAPHY型に変換することでメートルをそんまま扱えるのでテク ニックとして覚えておくと便利 SELECT id,ST_AsText(g),ST_Distance(g, ST_SetSRID(ST_MakePoint(134.97,34.999),6668)::geography) FROM t1 WHERE ST_Distance(g, ST_SetSRID(ST_MakePoint(134.97, 34.999),6668)::geography)<3000; id | | st_distance ----+----------------------------+--------------- 1 | POINT(135 35) | 2740.90786494 4 | POLYGON((135 35,135.0(略) | 2740.90786494 5 | POINT(134.95 35.02) | 2959.79779643

地理情報データの活用方法の例 交叉判定とその挙動の紹介 サンプルデータ SELECT id, ST_AsText(g) g FROM g1; id | g -------------------+------------------------------------------------------------------------------------------ basic_polygon | POLYGON((135 35,135.01 35,135.01 35.01,135 35.01,135 35)) line_cross1 | LINESTRING(134.995 35.009,135.003 35.007) line_cross2 | LINESTRING(134.995 35.003,135.012 35.007) line_nocross_in | LINESTRING(135.004 35.001,135.007 35.003) line_nocross_out | LINESTRING(134.998 35.011,135.008 35.012) polygon_inside | POLYGON((135.0085 35.008,135.009 35.008,135.009 35.009,135.0085 35.009,135.0085 35.008)) polygon_outside | POLYGON((135.011 35.008,135.015 35.008,135.015 35.012,135.011 35.012,135.011 35.008)) polygon_cross | POLYGON((134.997 34.998,135.003 34.998,135.003 35.002,134.997 35.002,134.997 34.998)) polygon_lineshare | POLYGON((135 35,135.01 35,135.01 34.996,135 34.996,135 35)) cross2 cross1 out in inside outside cross lineshare

地理情報データの活用方法の例 交叉判定関数の色々 ST_Crosses() ST_Intersects() ST_Contains() ST_Covers() ST_Disjoint() ST_Equals() ST_Overlaps() ST_Touches() ST_Within() SELECT a.id, b.id, ST_Crosses(a.g, b.g) crosses, ST_Intersects(a.g, b.g) intersects, ST_Contains(a.g, b.g) contains, ST_Covers(a.g, b.g) covers, ST_Disjoint(a.g, b.g) disjoint, ST_Equals(a.g, b.g) equals, ST_Overlaps(a.g, b.g) ovrlaps, ST_Touches(a.g, b.g) touches, ST_Within(a.g, b.g) witin FROM g1 a, g1 b WHERE a.id='basic_polygon'; 使用例

地理情報データの活用方法の例 • 交叉判定関数の挙動 SELECT a.id, b.id, ST_Crosses(a.g, b.g) crosses, ST_Intersects(a.g, b.g) intersects, ST_Contains(a.g, b.g) contains, ST_Covers(a.g, b.g) covers, ST_Disjoint(a.g, b.g) disjoint, ST_Equals(a.g, b.g) equals, ST_Overlaps(a.g, b.g) ovrlaps, ST_Touches(a.g, b.g) touches, ST_Within(a.g, b.g) witin FROM g1 a, g1 b WHERE a.id='basic_polygon'; id1 | id2 | crosses | intersects | contains | covers | disjoint | equals | ovrlaps | touches | witin ---------------+-------------------+---------+------------+----------+--------+----------+--------+---------+---------+------- basic_polygon | basic_polygon | f | t | t | t | f | t | f | f | t basic_polygon | line_cross1 | t | t | f | f | f | f | f | f | f basic_polygon | line_cross2 | t | t | f | f | f | f | f | f | f basic_polygon | line_nocross_in | f | t | t | t | f | f | f | f | f basic_polygon | line_nocross_out | f | f | f | f | t | f | f | f | f basic_polygon | polygon_inside | f | t | t | t | f | f | f | f | f basic_polygon | polygon_outside | f | f | f | f | t | f | f | f | f basic_polygon | polygon_cross | f | t | f | f | f | f | t | f | f basic_polygon | polygon_lineshare | f | t | f | f | f | f | f | t | f

地理情報データの活用：応用例 • メッシュデータ(5338-13-86-1～4)の例  このようなメッシュポリゴンが、メッシュIDと共に全国を埋め尽くして いる • ポリゴンデータの例(一部抜粋)  このようなポリゴンが全国に多数散らばっているデータ POLYGON ((138.4515632688455 35.487291658333326, 138.45156326888804 35.487312488319... POLYGON ((138.4518757689792 35.487916658333326, 138.45187576902188 35.487937488319... POLYGON ((138.44968827756946 35.48812498527777, 138.44968827761113 35.488145815263... POLYGON ((138.45218827512676 35.49145831638889, 138.45218827516965 35.491479156374... POLYGON ((138.4525007731528 35.491041657083336, 138.45250077319582 35.491062487069... POLYGON ((138.4525007727222 35.490833317222226, 138.4525007727653 35.4908541572083... POLYGON ((138.4528132711719 35.49062498777778, 138.4528132712151 35.49064581776389... POLYGON ((138.45281326901042 35.48958331847223, 138.45281326905365 35.489604158458... POLYGON ((138.45187577196876 35.48937498736111, 138.45187577201148 35.489395817347... POLYGON ((138.45000077208334 35.486041657083334, 138.450000772125 35.4860624870694,.. "533813861","POLYGON ((138.45624999999998 35.483333333333334, 138.45 35.4833333333... "533813862","POLYGON ((138.46249999999998 35.483333333333334, 138.45624999999998 3... "533813863","POLYGON ((138.45624999999998 35.4875, 138.45 35.4875, 138.45 35.49166... "533813864","POLYGON ((138.46249999999998 35.4875, 138.45624999999998 35.4875, 138...

地理情報データの活用：応用例 • メッシュ 5338-13-86-% と交叉するポリゴンを抽出する例 （ST_Intersects()を使用） SELECT m.key_code mesh_id, m.shape mesh, s.dataname ,s.g FROM chiri_mesh m LEFT OUTER JOIN shinsui_geom s ON (ST_Intersects(m.shape, s.g)) where m.key_code LIKE '53381386%';

地理情報データの活用：応用例 • メッシュ 5338-13-86-% と交叉するポリゴンを抽出し、占有割合を 算出する例(加えて、ST_Area()を使用) 実行結果 SELECT m.key_code mesh_id, SUM(ST_Area(ST_Intersection(m.shape, s.g))) / ST_Area(m.shape) wariai FROM chiri_mesh m LEFT OUTER JOIN shinsui_geom s ON (ST_Intersects(m.shape, s.g)) WHERE m.key_code LIKE '53381386%' GROUP BY m.key_code; mesh_id | wariai -----------+--------------------- 533813861 | 0.38508134563798374 533813862 | 0.5574807596776484 533813863 | 0.5151222089266544 533813864 | 0.7898798999943218

爆速DB Powered by PG-Strom • PG-Stromも、PostGIS関数に一部対応しています！ • 例えば、(複雑な形の)ポリゴン内に含まれる点の抽出など、 PG-Stromを使用することで数倍早く実行される実験結果が出 ています

爆速DB Powered by PG-Strom • PG-Stromが対応しているPostGIS関数 • st_makepoint() • st_point() • st_makepoint() • st_makepoint() • st_setsrid() • st_distance() • st_dwithin() • st_relate() • st_contains() • st_crosses() • st_linecrossingdirection() https://heterodb.github.io/pg-strom/ja/ref_devfuncs/#postgis ※まだまだ対応関数が少ないです。使う人が増えると対応してくれるかもしれないので、ぜひ皆さんお 試しください！

まとめ ぜひ データベースで 地理情報を 扱ってみましょう

One more thing… • 東大駒場にて「変わる高さ 動く大地 ～ 測 量に見せられた人々の物語」展 開催中 • 2024/9/28～11/24 ※本日行ってきました