GIS PostGIS 日付変更線マタギ問題

2024-12-032024
92.3% (337.8/366)
December
9.1% (2.8/31)
Week 49
26.1% (1.8/7)
Day 03 Tue
83.2% (19.9/24)

訪問者

total: 8403
today: 1
yesterday: 0
now: 1

↑

更新

MenuBar

何が問題?
検証
[方法1] 東西に二重に登録
[方法2] MULTIPOLYGON
[方法3] MULTIPOLYGON 簡易版
補遺データが POINT で、検索範囲が矩形の場合
補遺 POLYGON((...),(...)) について

何が問題? †

日付変更線をまたぐような矩形（①）を Postgis はうまく扱えない
ぢゃあ、ということで、無理やり左側の経度が右側の経度より小さくするような矩形(②、③)を作ってもうまくいかない

↑

検証 †

検証用のテーブルとデータ

CREATE TABLE geom_table (id bigserial, geom geometry(polygon, 4326));

insert into geom_table (geom) values (ST_GeomFromText('POLYGON((170 40, -170 40, -170 30, 170 30, 170 40))', 4326));
insert into geom_table (geom) values (ST_GeomFromText('POLYGON((170 40, 190 40, 190 30, 170 30, 170 40))', 4326));
insert into geom_table (geom) values (ST_GeomFromText('POLYGON((-190 40, -170 40, -170 30, -190 30, -190 40))', 4326));

# SELECT id, ST_AsText(geom) as area from geom_table;
 id |                        area
----+----------------------------------------------------
  1 | POLYGON((170 40,-170 40,-170 30,170 30,170 40))
  2 | POLYGON((170 40,190 40,190 30,170 30,170 40))
  3 | POLYGON((-190 40,-170 40,-170 30,-190 30,-190 40))
(3 rows)

どういうことが起きるか?

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(175 35)', 4326)) as distance from geom_table;
 id |                        area                        | distance
----+----------------------------------------------------+----------
  1 | POLYGON((170 40,-170 40,-170 30,170 30,170 40))    |        5
  2 | POLYGON((170 40,190 40,190 30,170 30,170 40))      |        0
  3 | POLYGON((-190 40,-170 40,-170 30,-190 30,-190 40)) |      345
(3 rows)

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(185 35)', 4326)) as distance from geom_table;
 id |                        area                        | distance
----+----------------------------------------------------+----------
  1 | POLYGON((170 40,-170 40,-170 30,170 30,170 40))    |       15
  2 | POLYGON((170 40,190 40,190 30,170 30,170 40))      |        0
  3 | POLYGON((-190 40,-170 40,-170 30,-190 30,-190 40)) |      355
(3 rows)

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(-185 35)', 4326)) as distance from geom_table;
 id |                        area                        | distance
----+----------------------------------------------------+----------
  1 | POLYGON((170 40,-170 40,-170 30,170 30,170 40))    |       15
  2 | POLYGON((170 40,190 40,190 30,170 30,170 40))      |      355
  3 | POLYGON((-190 40,-170 40,-170 30,-190 30,-190 40)) |        0
(3 rows)

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(-175 35)', 4326)) as distance from geom_table;
 id |                        area                        | distance
----+----------------------------------------------------+----------
  1 | POLYGON((170 40,-170 40,-170 30,170 30,170 40))    |        5
  2 | POLYGON((170 40,190 40,190 30,170 30,170 40))      |      345
  3 | POLYGON((-190 40,-170 40,-170 30,-190 30,-190 40)) |        0
(3 rows)

つまり
N35E175 N35E185 N35W175 N35W185
① ☓ ☓ ☓ ☓
② ○ ○ ☓ ☓
③ ☓ ☓ ○ ○
(☓: 矩形は該当座標を「含まない」と判定、○: 矩形は該当座標を「含む」と判定)
①の形式の「正しい」矩形の指定の仕方は、Postgis ではサポートされない !!

	N35E175	N35E185	N35W175	N35W185
①	☓	☓	☓	☓
②	○	○	☓	☓
③	☓	☓	○	○

↑

[方法1] 東西に二重に登録 †

日付変更線をまたぐような矩形(①)を Postgis に登録する場合には、②と③の矩形で二重に登録すればいい
東半球の POINT は②のみにヒット
西半球の POINT は③のみにヒット

日付変更線をまたぐような矩形と比較する場合には、一律東側に丸めて検索。②にヒットして③にはヒットしない

# SELECT id, ST_AsText(geom) as area,
  ST_Distance(geom, ST_GeomFromText('Polygon((175 35,185 35,185 30,175 30,175 35))', 4326)) as distance from geom_table;
 id |                        area                        | distance
----+----------------------------------------------------+----------
  1 | POLYGON((170 40,-170 40,-170 30,170 30,170 40))    |        5
  2 | POLYGON((170 40,190 40,190 30,170 30,170 40))      |        0
  3 | POLYGON((-190 40,-170 40,-170 30,-190 30,-190 40)) |      345
(3 rows)

↑

[方法2] MULTIPOLYGON †

[方法1] は、同じデータを二重に登録することになりキモイ

POLYGON ではなく MULTIPOLYGON を使って、矩形を日付変更線で分割してあげればデータを一意にできる

CREATE TABLE geom2_table (id bigserial, geom geometry(multipolygon, 4326));

insert into geom2_table (geom)
values (ST_GeomFromText('MULTIPOLYGON(
    ((170 40, 180 40, 180 30, 170 30, 170 40)),
    ((-180 40, -170 40, -170 30, -180 30, -180 40)))'
  , 4326));

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(-175 35)', 4326)) as distance from geom2_table;
 id |                                               area                                               | distance
----+--------------------------------------------------------------------------------------------------+----------
  1 | MULTIPOLYGON(((170 40,180 40,180 30,170 30,170 40)),((-180 40,-170 40,-170 30,-180 30,-180 40))) |        0
(1 row)

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(175 35)', 4326)) as distance from geom2_table;
 id |                                               area                                               | distance
----+--------------------------------------------------------------------------------------------------+----------
  1 | MULTIPOLYGON(((170 40,180 40,180 30,170 30,170 40)),((-180 40,-170 40,-170 30,-180 30,-180 40))) |        0
(1 row)

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(180 35)', 4326)) as distance from geom2_table;
 id |                                               area                                               | distance
----+--------------------------------------------------------------------------------------------------+----------
  1 | MULTIPOLYGON(((170 40,180 40,180 30,170 30,170 40)),((-180 40,-170 40,-170 30,-180 30,-180 40))) |        0
(1 row)

日付変更線上の点も矩形に含まれる(距離0)と認識される

↑

[方法3] MULTIPOLYGON 簡易版 †

[方法2] では、日付変更線で分割するのがメンドイ

① を ②と③ の MULTIPOLYGON にしてもよい

insert into geom2_table (geom)
values (ST_GeomFromText('MULTIPOLYGON(
    ((170 40, 190 40, 190 30, 170 30, 170 40)),
    ((-190 40, -170 40, -170 30, -190 30, -190 40))
  , 4326));

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(-175 35)', 4326)) as distance from geom2_table;
 id |                                               area                                               | distance
----+--------------------------------------------------------------------------------------------------+----------
  1 | MULTIPOLYGON(((170 40,180 40,180 30,170 30,170 40)),((-180 40,-170 40,-170 30,-180 30,-180 40))) |        0
  2 | MULTIPOLYGON(((170 40,190 40,190 30,170 30,170 40)),((-190 40,-170 40,-170 30,-190 30,-190 40))) |        0
(2 rows)

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(175 35)', 4326)) as distance from geom2_table;
 id |                                               area                                               | distance
----+--------------------------------------------------------------------------------------------------+----------
  1 | MULTIPOLYGON(((170 40,180 40,180 30,170 30,170 40)),((-180 40,-170 40,-170 30,-180 30,-180 40))) |        0
  2 | MULTIPOLYGON(((170 40,190 40,190 30,170 30,170 40)),((-190 40,-170 40,-170 30,-190 30,-190 40))) |        0
(2 rows)

SELECT id, ST_AsText(geom) as area, ST_Distance(geom, ST_GeomFromText('Point(180 35)', 4326)) as distance from geom2_table;
 id |                                               area                                               | distance
----+--------------------------------------------------------------------------------------------------+----------
  1 | MULTIPOLYGON(((170 40,180 40,180 30,170 30,170 40)),((-180 40,-170 40,-170 30,-180 30,-180 40))) |        0
  2 | MULTIPOLYGON(((170 40,190 40,190 30,170 30,170 40)),((-190 40,-170 40,-170 30,-190 30,-190 40))) |        0
(2 rows)

きちんと分割しても、いい加減に分割しても同じ

↑

補遺データが POINT で、検索範囲が矩形の場合 †

データは、日付変更線をまたぐことはない (POINTだから)
したがって、検索範囲を日付変更線で分割して MULTIPOLYGON にしてやる必要がある。
まとめると
- データが矩形、検索範囲がPOINT の場合
  - → データが日付変更線をまたぐときには、２重に登録
- データが矩形、検索範囲が矩形の場合
  - → データが日付変更線をまたぐときには、２重に登録
  - → 検索範囲が日付変更線をまたぐときには、検索範囲を東側に丸めた POLYGON で検索
- データがPOINT、検索範囲が矩形の場合
  - → 検索範囲を日付変更線で分割して MULTIPOLYGON にして検索
どの場合でも一律、日付変更線で分割した MULTIPOLYGON を使えば、東経185度などの変則的な変換をしなくてもいい。遅いかもしれんけど

↑

補遺 POLYGON((...),(...)) について †

POLYGON((...),(...)) は、MULTIPOLYGON ではなく、中空のドーナッツ型を示す。

insert into geom_table (geom) values (
ST_GeomFromText('POLYGON((-20 20, 20 20, 20 -20, -20 -20, -20 20),(-10 10, 10 10, 10 -10, -10 -10, -10 10))', 4326));

SELECT id, ST_AsText(geom) as area from geom_table where ST_Distance(geom, ST_GeomFromText('Point(-15 0)', 4326)) = 0;
 id |                                        area
----+------------------------------------------------------------------------------------
  4 | POLYGON((-20 20,20 20,20 -20,-20 -20,-20 20),(-10 10,10 10,10 -10,-10 -10,-10 10))
(1 row)

SELECT id, ST_AsText(geom) as area from geom_table where ST_Distance(geom, ST_GeomFromText('Point(15 0)', 4326)) = 0;
 id |                                        area
----+------------------------------------------------------------------------------------
  4 | POLYGON((-20 20,20 20,20 -20,-20 -20,-20 20),(-10 10,10 10,10 -10,-10 -10,-10 10))
(1 row)

SELECT id, ST_AsText(geom) as area from geom_table where ST_Distance(geom, ST_GeomFromText('Point(0 0)', 4326)) = 0;
 id | area
----+------
(0 rows)

GIS

GIS PostGIS 日付変更線マタギ問題

目次

訪問者

更新

最新の10件

人気の10件

MenuBar

何が問題? †

検証 †

[方法1] 東西に二重に登録 †

[方法2] MULTIPOLYGON †

[方法3] MULTIPOLYGON 簡易版 †

補遺データが POINT で、検索範囲が矩形の場合 †

補遺 POLYGON((...),(...)) について †

GIS PostGIS 日付変更線マタギ問題

目次

訪問者

更新

最新の10件

人気の10件

MenuBar

何が問題? †

検証 †

[方法1] 東西に二重に登録 †

[方法2] MULTIPOLYGON †

[方法3] MULTIPOLYGON 簡易版 †

補遺 データが POINT で、検索範囲が矩形の場合 †

補遺 POLYGON((...),(...)) について †

補遺データが POINT で、検索範囲が矩形の場合 †