18. ジオグラフィ¶
座標が「地理座標系」または「緯度/経度」のデータを持つことは極めて一般的です。
地理座標系は、メルカトル、UTM、州平面、座標と違い、デカルト座標ではありません。地理座標系では、平面上にプロットするような原点からの線形距離を表現しない座標系です。これらの**球面座標**は、地球上の角度の座標を表します。球面座標では、ポイントは基準子午線からの回転角度 (経度) と赤道からの角度 (緯度) で指定されます。
地理座標をデカルト座標系の近似のように扱えますし、空間計算を続行できます。しかしながら、距離、長さ、面積の計測は**無意味**です。球面座標系は**角度**の距離を計測し、単位は「度」となるためです。さらに、インデックスとインタセクトや包含といった真偽テストからの近似的な結果が、ひどく悪くなる可能性があります。ポイント間の距離は、極 (北極や南極) や日付変更線といった問題の領域では大きくなります。
たとえば、ロサンゼルスとパリの座標です。
ロサンゼルス:
POINT(-118.4079 33.9434)パリ:
POINT(2.3490 48.8533)
通常のPostGISデカルト距離の計算である ST_Distance(geometry, geometry) を使ってロサンゼルスとパリの間の距離の計算は次のように行います。SRIDが4326で地理空間参照系を宣言していることに注意して下さい。
SELECT ST_Distance(
'SRID=4326;POINT(-118.4079 33.9434)'::geometry, -- Los Angeles (LAX)
'SRID=4326;POINT(2.5559 49.0083)'::geometry -- Paris (CDG)
);
121.898285970107
わっ! 122! でもこれはどういう意味でしょう?
空間参照系 4326 の単位は度です。結果は 122 度でした。でも (二度目)、これはどういう意味でしょう?
球面では、「1平方度」は非常に可変的です。赤道から離れるにつれ小さくなります。極に向かって進むにつれ、地球上の子午線 (垂直線) が他の子午線と近くなっていくと考えて下さい。122度という距離は、全く*意味がありません*。無意味な数字なのです。
意味のある距離を計算するためには、地理座標を近似的なデカルト座標として扱わずに、本当の球面座標系として扱わなければなりません。ポイント間の距離について、本当の球面上の経路である大円の一部に沿った計測が必要です。
PostGISはこの機能を「ジオグラフィ (geography)」型を通して提供しています。
注釈
空間データベースが異なると、「ジオグラフィの処理」のアプローチも異なります
ORACLEは、SRIDが地理座標系の時に、透過的な地理座標系計算によって、違いを覆い隠そうとします。
SQL Serverでは、デカルト座標では"STGeometry"、地理座標では"STGeography"を、それおれ使います。
Informix Spatial は、純粋なデカルト座標のInformix拡張で、また、Informix Gedeticが純粋な地理座標拡張です。
PostGISは「ジオメトリ」と「ジオグラフィ」の二つの型を使っていて、SQL Serverと近いです。
では、geometry のかわりに geography を使って、ロサンゼルスとパリの間の距離の計測を、もう一度やってみます。
SELECT ST_Distance(
'SRID=4326;POINT(-118.4079 33.9434)'::geography, -- Los Angeles (LAX)
'SRID=4326;POINT(2.5559 49.0083)'::geography -- Paris (CDG)
);
9124665.27317673
大きい値ですね! 「ジオグラフィ」計算からの返り値は全て**メートル**単位で、結果が 9125km です。
古いバージョンのPostGISでは、ST_Distance_Spheroid(point, point, measurement) を使った球面上の非常に基本的な計算がサポートされていました。しかしながら、ST_Distance_Spheroid は実質的に制限されます。この関数はポイント同士でのみ有効であり、極や日付変更線を越えるインデックスのサポートはありません。
「ロサンゼルスからパリまでのフライトでどれだけアイスランドに近づくか」等の問題が出る時に、ポイント以外のジオメトリのサポートの必要性が非常にはっきりしました。
デカルト平面 (紫線) で地理座標を使うと、確かに、*非常に*誤った答えが得られます! 大円ルート (赤線) を使うと正しい答えが得られます。「ジオグラフィ」を使ってLAX-CDG (ロサンゼルス-ド・ゴール) のフライトをラインストリングに変換して、アイスランドのあるポイントまでの距離を計算すると、メートル単位の正しい答が得られます。
SELECT ST_Distance(
ST_GeographyFromText('LINESTRING(-118.4079 33.9434, 2.5559 49.0083)'), -- LAX-CDG
ST_GeographyFromText('POINT(-22.6056 63.9850)') -- Iceland (KEF)
);
502454.906643729
LAX-CDGルートでアイスランドに最も接近する距離 (アイスランドの国際空港からの距離とします) は比較的小さな502kmです。
ジオグラフィ座標の処理へのデカルトのアプローチは、日付変更線をクロスする地物について、完全に分解します。ロサンゼルスから東京までの最短大円ルートは太平洋とクロスします。デカルトでの最短ルートは大西洋とインド洋とクロスします。
SELECT ST_Distance(
ST_GeometryFromText('Point(-118.4079 33.9434)'), -- LAX
ST_GeometryFromText('Point(139.733 35.567)')) -- NRT (Tokyo/Narita)
AS geometry_distance,
ST_Distance(
ST_GeographyFromText('Point(-118.4079 33.9434)'), -- LAX
ST_GeographyFromText('Point(139.733 35.567)')) -- NRT (Tokyo/Narita)
AS geography_distance;
geometry_distance | geography_distance
-------------------+--------------------
258.146005837336 | 8833954.76996256
18.1. ジオグラフィを使う¶
ジオメトリデータをジオグラフィテーブルにロードするために、ジオメトリは最初に EPSG:4326 (経度/緯度) に投影変換する必要があります。それから、ジオグラフィに変換します。 ST_Transform(geometry,srid) 関数は座標で地理座標系に変換し、 Geography(geometry) 関数または ::geography 後置詞でジオグラフィに「キャスト」します。
CREATE TABLE nyc_subway_stations_geog AS
SELECT
ST_Transform(geom,4326)::geography AS geog,
name,
routes
FROM nyc_subway_stations;
ジオグラフィテーブルで空間インデックスを構築することは、全くジオメトリと同じです:
CREATE INDEX nyc_subway_stations_geog_gix
ON nyc_subway_stations_geog USING GIST (geog);
違いは極または日付変更線をまたぐ場合にでてきます。ジオグラフィインデックスは極または日付変更線をまたぐクエリを正しく処理しますが、ジオメトリは正しくは処理しません。
次に、エンパイアステートビルから500メートル圏内にある全ての地下鉄の駅を検索するクエリを示します。
ジオグラフィ型用のネイティブ関数だけが少数存在します:
ST_AsText(geography) は
text型を返しますST_GeographyFromText(text) は
geography型を返しますST_AsBinary(geography) は
bytea型を返しますST_GeogFromWKB(bytea) は
geography型を返しますST_AsSVG(geography) は
text型を返しますST_AsGML(geography) は
text型を返しますST_AsKML(geography) は
text型を返しますST_AsGeoJson(geography) は
text型を返しますST_Distance(geography, geography) は
double型を返しますST_DWithin(geography, geography, float8) は
boolean型を返しますST_Area(geography) は
double型を返しますST_Length(geography) は
double型を返しますST_Covers(geography, geography) は
boolean型を返しますST_CoveredBy(geography, geography) は
boolean型を返しますST_Intersects(geography, geography) は
boolean型を返しますST_Buffer(geography, float8) は
geography1 型を返しますST_Intersection(geography, geography) は
geography1 型を返します
18.2. ジオグラフィテーブルの生成¶
ジオグラフィカラムを持つテーブルの新規作成のSQLはジオメトリテーブルの生成とよく似ています。しかしながら、ジオグラフィはテーブル生成時にオブジェクトタイプを直接指定する機能が含まれます。例を示します:
CREATE TABLE airports (
code VARCHAR(3),
geog GEOGRAPHY(Point)
);
INSERT INTO airports
VALUES ('LAX', 'POINT(-118.4079 33.9434)');
INSERT INTO airports
VALUES ('CDG', 'POINT(2.5559 49.0083)');
INSERT INTO airports
VALUES ('KEF', 'POINT(-22.6056 63.9850)');
テーブル定義では GEOGRAPHY(Point) は空港データタイプをポイントとして指定します。新しいジオグラフィフィールドは geometry_columns ビューには登録されません。かわりに geography_columns と呼ばれるビューに登録されます。
SELECT * FROM geography_columns;
f_table_name | f_geography_column | srid | type
--------------------------+--------------------+------+----------
nyc_subway_stations_geog | geog | 0 | Geometry
airports | geog | 4326 | Point
注釈
上の出力から一部のカラムが省略されました。
18.3. ジオメトリへのキャスト¶
ジオグラフィ型の基本関数が多数のユースケースを処理できる一方、ジオメトリ型にしか対応しない関数にアクセスする必要がある場合があります。幸運なことに、ジオメトリとジオグラフィは相互に変換できます。
PostgreSQLのキャストに関する構文規則では、キャストしたい値の後ろに ::型名 を追加します。よって 2::text は 2 という数値を '2' という文字列に変換します。また、'POINT(0 0)'::geometry はポイントの文字列表現をジオメトリのポイントに変換します。
ST_X(point) 関数はジオメトリ型だけに対応しています。では、ジオグラフィからX値を読むにはどうすればいいでしょうか?
SELECT code, ST_X(geog::geometry) AS longitude FROM airports;
code | longitude
------+-----------
LAX | -118.4079
CDG | 2.5559
KEF | -21.8628
By appending ::geometry to our geography value, we convert the object to a geometry with an SRID of 4326. From there we can use as many geometry functions as strike our fancy. But, remember -- now that our object is a geometry, the coordinates will be interpreted as Cartesian coordinates, not spherical ones.
18.4. なぜジオグラフィを使うのか/使わないのか¶
地理座標は普遍的に受け入れられる座標で、誰もが緯度/経度がどういう意味かを理解していますが、UTM座標がどういう意味かを理解している人はまずいません。なのになぜ地理座標系を常に使うことにならないのでしょうか?
まず、前述の通り、ジオグラフィ型を直接サポートする関数は (今のところ) ごく少数です。ジオグラフィ型の制限の回避に多くの時間を費やす場合もあります。
それと、球面上の計算は、デカルト計算よりも非常に計算量が多くなります。たとえば、デカルト空間の距離に関する式 (ピタゴラスの定理) は、sqrt() を1回呼びます。球面上の距離に関する式 (半正矢関数) は、sqrt() を2回呼び、arctan() を1回呼び、sin() を4回呼び、cos() を2回呼びます。三角関数は非常に計算コストがかかり、球面計算には多数の三角関数が使われます。
結論は?
データが地理的にコンパクトな場合 (州、郡、市に含まれる場合)、ジオメトリ型とデカルト投影を使います。http://epsg.io/ サイトを見て、使用可能な空間参照系を選択するために、お住いの地域の名前を入力してみて下さい。
地理的に分散したデータセットで距離を測定する必要がある場合には (世界の大部分をカバーする場合)、ジオグラフィ型を使います。「ジオグラフィ」の作業でアプリケーションの複雑性が節約できるので、あらゆるパフォーマンス上の問題を埋め合わせてくれます。そして、「ジオメトリ」は多くの機能制限を埋め合わす能力があります。
18.5. 関数リスト¶
ST_Distance(geometry, geometry): ジオメトリ型については、二つのジオメトリの2次元デカルト距離の最小値 (空間参照系に基づく) を空間参照系の単位で返します。ジオグラフィ型については、デフォルトでは二つのジオメトリの球面の最小距離をメートル単位で返します。
ST_GeographyFromText(text): Well-Known Text (WKT) 表現または拡張 WKT (EWKT) 表現から、ジオグラフィを返します。
ST_Transform(geometry, srid): 整数パラメータで指定される SRID に座標変換した新しいジオメトリを返します。
ST_X(point): ポイントのX座標値を返します。有効でないなら NULL を返します。入力はポイントでなければなりません。
ST_Azimuth(geography_A, geography_B): AからBの方向をラジアン単位で返します。
ST_DWithin(geography_A, geography_B, R): AがBのRメートル以内であるなら TRUE を返します。
脚注
