LoginSignup
1
0

More than 1 year has passed since last update.

@ak-sakatoku(Akira Sakatoku)inKDDI株式会社

SnowflakeがSRIDに対応したので勉強した

Last updated at Posted at 2023-05-02

「つい最近のアップデートでSnowflakeのGEOMETRY型がSRIDに対応したぞ！」と周りのデータエンジニアさんが盛り上がっていたんですが…実はSRIDをよく知らなかったので勉強してみました。
Snowflake Documentation | What's New - April 2023

SRIDとは

Wikipediaによると、

A spatial reference system (SRS) or coordinate reference system (CRS) is a framework used to precisely measure locations on the surface of the Earth as coordinates.
A Spatial Reference System Identifier (SRID) is a unique value used to unambiguously identify projected, unprojected, and local spatial coordinate system definitions. These coordinate systems form the heart of all GIS applications.

とのことで、これを日本語にすると

空間参照系(SRS)もしくは座標参照系(CRS)は地球の表面上の座標を正確に測定するために使われるフレームワークです。
空間参照系IDは投影座標系、非投影座標系、およびローカル座標系の定義を識別するために使われるユニークな値です。これらの座標系はあらゆるGISアプリケーションの中核です。

という感じでしょうか。初心者には用語一つ一つが難しいんですが、

  • 投影座標系：三次元の座標を地図上の二次元の座標に変換した座標系
  • 非投影座標系：投影していない、つまり三次元の座標そのままの座標系
  • ローカル座標系：ある場所を基準として地球規模ではない限られた範囲を対象にした座標系

ってことかな？

Wikipediaの例では、それぞれのSRSの定義に使われるパラメータがまちまちですが、例えば

  • 測地系(Datum)：地球をどんな形で表してどんな座標系を重ねるか。GRS80楕円体＋ITRF94座標系、WGS84楕円体＋WGS84座標系、…
  • 基準点：グリニッジ天文台を基準にする、みたいな話
  • 測量単位(Unit of Measure)：度、メートル、…

などをパラメータとして定義しているようです。

どんなSRSがあるか調べてみる

SRIDは正式には出典元を示す文字列(オーソリティ)と整数のIDを並べてEPSG:4326みたいに記載するそうです。この例ではEPSG1が出典元。
ただ、EPSGが管理しているものがデファクトスタンダードになっているので、オーソリティを省略して単に4326と記載することも多いようです。
とはいえ、EPSGが管理していない野良のSRIDもあるとのこと。

EPSGが管理しているSRIDのデータベース、SRIDレジストリは以下のサイト。
EPSG Geodetic Parameter Dataset
どうでもいいけど検索に癖があるなあ…

4326で検索してみると、
WGS 84 (EPSG:4326)
SRIDはそのまんまEPSG:4326。WGS84ってどこかで見たことある！

TokyoとRectangularで検索してみると、
JGD2000 / Japan Plane Rectangular CS IX (EPSG:2451)
SRIDはEPSG:2451。WGS84と色々違うけど測量単位がメートル！

Africa、west、Capeで検索してみると、
Cape / Lo17 (EPSG:22277)
SRIDはEPSG:22277。南アフリカの測量で使う…みたいなことが書かれていますが、軸の向きが上記2つと違いますね。

135, 35という値を渡されたら日本のどこかを指す経度と緯度？と思ってしまいますが、その値がどのSRSにおけるものなのか次第で差している場所が変わりそうということが分かりました。

SnowflakeでSRIDに触れてみる

なんとなくどういうものか分かったところで、とりあえずリリースノート通りのSQL文を実行してみます。

-- X-Smallウェアハウスを選択
USE WAREHOUSE COMPUTE_WH;

-- テスト用のデータベースとスキーマを作成する
CREATE OR REPLACE DATABASE DEV_TEST;
USE DATABASE DEV_TEST;
CREATE OR REPLACE SCHEMA TEST_SRID;
USE SCHEMA TEST_SRID;

-- リリースノート通りのSQL文を発行してみる
SELECT TO_GEOMETRY('POINT(1820.12 890.56)', 4326);
-- { "coordinates": [ 1.820120000000000e+03, 8.905599999999999e+02 ], "type": "Point" }

-- デフォルトの設定のままだとSRIDを変えても出力は一文字も変わらない
SELECT TO_GEOMETRY('POINT(1820.12 890.56)', 3097);
-- { "coordinates": [ 1.820120000000000e+03, 8.905599999999999e+02 ], "type": "Point" }

-- 出力フォーマットをEWKTに変換してみる
ALTER SESSION SET GEOMETRY_OUTPUT_FORMAT = EWKT;
SELECT TO_GEOMETRY('POINT(1820.12 890.56)', 4326);
SELECT TO_GEOMETRY('POINT(1820.12 890.56)', 3097);
-- SRID=4326;POINT(1820.12 890.56)
-- SRID=3097;POINT(1820.12 890.56)

出力フォーマットをEWKT2にしてみると、値にSRIDが添えられていることが確認できました！

地理空間関数を使ってみる

GEOMETRY型の値の書き方が分かったので、何か地理空間関数を使ってみましょう。ST_DISTANCEを使えば点と点の距離を測れそう。

-- ST_DISTANCEを試してみる
SELECT ST_DISTANCE(
  TO_GEOMETRY('POINT(135 38)', 4326),
  TO_GEOMETRY('POINT(136 39)', 4326)
) AS DISTANCE;
-- DISTANCE
-- 1.414213562

ふむふむ。SRIDを書き換えて、SRSが違うデータをST_DISTANCEに渡してあげるとどうなるのかな？

-- ST_DISTANCEを試してみる
SELECT ST_DISTANCE(
  TO_GEOMETRY('POINT(135 38)', 4326),
  TO_GEOMETRY('POINT(135 38)', 4612)
) AS DISTANCE_BETWEEN_SRS;
-- 100388 (P0000): Incompatible SRID: 4326 and 3097

なるほど、地理空間関数はSRSが違うデータの比較・計算に対応していないからエラーが出るみたいですね。

GEOMETRY型を含むテーブルを作る

GEOMETRY型を含むテーブルを作ってテストデータを格納してみましょう。テーブルに格納するときは例えばEWKTで書いたものをそのまま投げ込めばいいみたいです。

-- テストデータを格納するテーブルを作成する
CREATE OR REPLACE TABLE TEST_GEO (
id INTEGER,
geo GEOMETRY
);
-- テストデータを手動で格納していく。最後の1行だけSRIDを変えてみる
INSERT INTO TEST_GEO(id, g)
VALUES
(1, 'SRID=4326;POINT(134 38)'),
(2, 'SRID=4326;POINT(135 38)'),
(3, 'SRID=4326;POINT(136 38)'),
(4, 'SRID=4326;POINT(134 39)'),
(5, 'SRID=4326;POINT(135 39)'),
(6, 'SRID=4326;POINT(136 39)'),
(7, 'SRID=4326;POINT(134 40)'),
(8, 'SRID=4326;POINT(135 40)'),
(9, 'SRID=4326;POINT(136 40)'),
(10, 'SRID=4326;POINT(134 38)'),
(11, 'SRID=4612;POINT(134 38)');
-- テストデータが格納されたことを確認する
SELECT * FROM TEST_GEO;

テーブルに格納したデータのSRIDを確認する

ST_SRIDを使えばSRIDを取得できるらしいので、これで先ほどのテーブルに格納したデータのSRIDを確認してみましょう。

-- SRIDだけを確認してみる
SELECT DISTINCT ST_SRID(g) AS srid FROM TEST_GEO;
-- SRID
-- 4326
-- 4612

テーブルに格納したデータに地理空間関数を使ってみる

とりあえず点と点の距離を測ってみましょう。

SELECT a.id, ST_DISTANCE(a.g, b.g) AS distance
FROM TEST_GEO a
JOIN (
  -- 1行目を基準点として計算する
  SELECT g
  FROM TEST_GEO
  WHERE id = 1
) b;
-- 100388 (P0000): Incompatible SRID: 4612 and 4326

むむ。SRSが違うデータが混ざっているからエラーが。先ほどのST_SRIDを使ってSRSが違うデータを除外してみましょう。

SELECT a.id, ST_DISTANCE(a.g, b.g) AS distance
FROM TEST_GEO a
JOIN (
  -- 1行目を基準点として計算する
  SELECT g
  FROM TEST_GEO
  WHERE id = 1
) b
ON ST_SRID(a.g) = ST_SRID(b.g);
-- ID	DISTANCE
-- 1	0
-- 2	1
-- 3	2
-- 4	1
-- 5	1.414213562
-- 6	2.236067977
-- 7	2
-- 8	2.236067977
-- 9	2.828427125
-- 10	0

無事に成功して点と点の距離を計算できました！このケースでは測量単位がメートルじゃないのに単純計算している、ので実用上はあまり意味ないかも。
あと、ドキュメントには以下の記載がありました。同じ列にSRSが違うデータが混ざっているとパフォーマンス最適化に制限が生じるようです。いちいちWHERE句で除外するのも面倒なので、できるだけ統一しておきたいですね。

GEOMETRY列には、異なるSRIDsを持つオブジェクトを挿入できます。列に複数のSRIDが含まれている場合は、重要なパフォーマンス最適化の一部が適用されません。これにより、特に地理空間述語で結合する場合に、クエリが遅くなる可能性があります。

SRIDを設定する

SRIDを手動で設定したくなる場合もありそうですよね。
SRIDがデータに明記されていないけど仕様書で分かっている場合、誤差が生じるのを承知でSRSを統一したい場合など。
ST_SETSRIDが使えそうです。

SELECT a.id, ST_DISTANCE(ST_SETSRID(a.g, b.refer_srid), b.g) AS distance, b.refer_srid
FROM TEST_GEO a
JOIN (
  -- 1行目を基準点として計算する
  SELECT g, ST_SRID(g) as refer_srid
  FROM TEST_GEO
  WHERE id = 1
) b;
-- ID	DISTANCE	REFER_SRID
-- 1	0			4326
-- 2	1			4326
-- 3	2			4326
-- 4	1			4326
-- 5	1.414213562	4326
-- 6	2.236067977	4326
-- 7	2			4326
-- 8	2.236067977	4326
-- 9	2.828427125	4326
-- 10	0			4326
-- 11	0			4326

この場合はすべての点のSRIDを4326として計算していて、結果的に11番目の点のSRIDが4612から変更されています。SRID：4326はWGS84、SRID：4612は日本測地系2000のことで、誤差はあるけどそこまで大きくない…らしいので、大きな問題はなさそう。
値を変換している訳ではなくて値に添えられているSRIDを変更しているだけなので、測量単位が度(緯度・経度)のSRSとメートルのSRSが混ざったデータでこういう計算をするととんでもないことになりそう！

SRIDを設定したデータを保存する

SRIDを設定したデータを何回も使うならテーブルに格納したくなりますね。
CTASクエリ3でやっちゃうのが楽そうです。

-- CTASクエリでSRIDを設定したデータをテーブルに保存する
CREATE OR REPLACE TABLE TEST_GEO2 (id, g)
AS SELECT id, ST_SETSRID(g, 4326) FROM TEST_GEO;
SELECT * FROM TEST_GEO2;

他の地理空間関数も使ってみよう

ST_WITHINを使うと、例えばある点がポリゴンの内側にあるかどうかを調べられるようです。
ある都道府県のデータだけフィルタするときに使えそう…試してみましょう。

-- ST_WITHINを試してみる
SELECT a.id, a.g
FROM TEST_GEO a
JOIN (
  SELECT
  TO_GEOMETRY('SRID=4326;POLYGON((135.976 39.569, 135.627 39.955, 135.107 39.924, 134.808 39.497, 134.955 38.998, 135.437 38.801, 135.891 39.055, 135.976 39.569))') AS area
) AS b
WHERE ST_SRID(a.g) = ST_SRID(b.area) AND ST_WITHIN(a.g, b.area);
-- ID	G
-- 5	SRID=4326;POINT(135 39)

ここでもSRSが違うとエラーが発生します。今回のテストデータはSRSが違うデータが混ざっているのでST_SRIDでSRIDが同じことを確認しています。
ST_WITHINの比較に用いているポリゴンは135.4, 39.4を中心にした半径0.6の正多角形で、その内側に1件のテストデータが見つかりました。

-- ST_WITHINを試してみる
SELECT a.id, a.g
FROM TEST_GEO a
JOIN (
  SELECT
  TO_GEOMETRY('SRID=4326;POLYGON((135.784 39.513, 135.551 39.77, 135.205 39.749, 135.005 39.465, 135.103 39.132, 135.424 39.001, 135.727 39.17, 135.784 39.513))') AS area
) AS b
WHERE ST_SRID(a.g) = ST_SRID(b.area) AND ST_WITHIN(a.g, b.area);
-- Query produced no results

ST_WITHINの比較に用いるポリゴンを半径0.4に変更すると、ちゃんと何も表示されなくなりました。

感想

SRIDはなんとなく分かりましたが、地理空間情報は定義や理論のレベルで専門性が高いなあというのが体感できました。
でもSnowflakeはサクッと触れて公式ドキュメントも充実しているので、GEOMETRY型のデータを試してみる、程度であれば全然困りませんでした。ありがとうSnowflake！
Snowflakeドキュメント | 地理空間データ型
Snowflakeドキュメント | 地理空間関数

  1. European Petroleum Survey Group(欧州石油調査グループ)。ヨーロッパの石油産業に由来する団体で2005年まで活動していた。石油を発掘するためには測地・測量が大事らしいです。 ↩

  2. Extended Well-Known Text。Wikipediaによると、PostGISの独自フォーマットだそうです。 ↩

  3. CREATE TABLE ... AS SELECT ...というクエリ。SELECTした結果で新しいテーブルを作成する。 ↩

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0