緯度経度の数字で絞るのと緯度経度の範囲(POLYGON)で絞るのとは違う

先日の shunyasu さんのトークをきっかけに、以下の日記を書きました。

sakaik.hateblo.jp

この中で、私が

> GeoHashも結局緯度経度で絞り込んでいることを考えると、すべてのクエリで同じ件数が得られるべきと考えるのが自然

と書いたのですが、これは誤りでした。
その後の shunyasuさんとの Twitter(X)でのやりとりの中で発覚しました。

両者の違い

　単純に緯度経度の数字の範囲に収まっているものを抽出するのと、その緯度経度の各点を頂点としたポリゴンの中に含まれるものを抽出するのは、一見すると同じことのように思えます。
ポリゴンだって、その緯度と経度の範囲を検索するわけですからね。
しかし、ポリゴンのある2つの頂点を結ぶ点がどこを通るのかを考えてみると、実は違うことに気づきます。

　学生の頃などに、メルカトル図法の地図を見ながら、「東京とサンフランシスコ（あるいはNYだったりパリだったり）の飛行機は、こんなに "遠回り"をしている（実はそれが最短距離）」という話を聞いたことがあると思います。今回の話題は、まさにそれでした。　もう一度書くのも面倒なので、以下のツイートを参照してください。

> GeoHashも結局緯度経度で絞り込んでいることを考えると、すべてのクエリで同じ件数が得られるべきと考えるのが自然

Geohashと緯度経度の結果は同じでした！たまにSRID 4326とズレる感じでしたね。ごく稀にSRID 0ともズレることがあり、これは完全に緯度経度を保存できてないからですかね...
— shunyasu (@ssuzuki67) 2025年12月23日

やっぱり、緯度経度の数字だけで範囲を絞る(A)のと緯度経度の4点を指定して地理座標上でその中に含まれる点を調べるのは違う操作です。
よく言う「日本から西海岸への最短距離」を例にすると、緯度経度数値内の範囲を調べる辺は(A)、地球上での2点を結ぶ辺は(B)なので、結果が変わるのも然りと。 pic.twitter.com/fE4SQ3hkeI
— 坂井恵(SAKAI Kei) (@sakaik) 2025年12月29日

別の部分を切り出してみると、こういう図でも説明可能かと思います。右図の赤い点が、緯度経度の数字だけで見ると含まれたり含まれなかったりしますね。
寧ろ、緯度経度の数字だけで見るケースと、デカルト（SRID=0）では(正しく処理すれば)一致するはず、という気がします。 pic.twitter.com/SSdLGclZY0
— 坂井恵(SAKAI Kei) (@sakaik) 2025年12月29日

ということで試してみた

　下図のような、(35 135)－(36 136) の四角形のポリゴンを考え、検索手法によりPOINTが内部に含まれるかどうかを判定したいと思います。

　ここでは下辺（35N で 135E から 136E を結ぶ辺）に注目して実験することにします。北緯35度を跨ぐ点をいくつか登録します。

CREATE TABLE points (p POINT SRID 4326);

INSERT INTO points VALUES (ST_GeomFromText('POINT(35.005 135.5)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(35.004 135.5)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(35.003 135.5)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(35.002 135.5)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(35.001 135.5)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(35 135.5)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(34.999999 135.5)',4326));

こんなデータが登録されます。

mysql> SELECT ST_AsText(p) FROM points;
+------------------------+
| ST_AsText(p)           |
+------------------------+
| POINT(35.005 135.5)    |
| POINT(35.004 135.5)    |
| POINT(35.003 135.5)    |
| POINT(35.002 135.5)    |
| POINT(35.001 135.5)    |
| POINT(35 135.5)        |
| POINT(34.999999 135.5) |
+------------------------+
7 rows in set (0.000 sec)

緯度経度の数値で絞り込む

　まず、単純に緯度経度の数値の大小で絞り込むクエリを確認してみましょう。

mysql> SELECT ST_AsText(p) FROM points
    ->  WHERE ST_Latitude(p) BETWEEN 35 AND 36
    ->    AND ST_Longitude(p) BETWEEN 135 AND 136;
+---------------------+
| ST_AsText(p)        |
+---------------------+
| POINT(35.005 135.5) |
| POINT(35.004 135.5) |
| POINT(35.003 135.5) |
| POINT(35.002 135.5) |
| POINT(35.001 135.5) |
| POINT(35 135.5)     |
+---------------------+
6 rows in set (0.000 sec)

　期待通り、北緯35度線よりも北にある点がすべて選択されました（"BETWEEN" なので境界線上の点も含まれています）。

ポリゴンによる絞り込み

　次に同じ範囲（のつもり）のポリゴンを作成して絞り込んでみます。

mysql> SELECT ST_AsText(p) FROM points 
    -> WHERE ST_WITHIN(p, ST_GeomFromText('POLYGON((35 135, 36 135, 36 136, 35 136, 35 135))',4326));
+---------------------+
| ST_AsText(p)        |
+---------------------+
| POINT(35.005 135.5) |
| POINT(35.004 135.5) |
| POINT(35.003 135.5) |
| POINT(35.002 135.5) |
+---------------------+
4 rows in set (0.000 sec)

　おや！？（←わざとらしい）　北緯35度、北緯35.001度の点が含まれなくなってしまいました。
これは地球上で 135Eと136Eを結ぶ「直線」は、その中間点(135.5E)に於いて 35.001Nと35.002Nの間を通っていることを意味します。
この狭い範囲でさえ、見事に西海岸行きの飛行機を再現できました。

まとめ

　メルカトル図法に飼い慣らされてしまった我々にとって、このST_WITHIN()の結果は直感に反するかもしれません。これが地理座標における「正しい範囲」と言えます。
一方、緯度経度の範囲で切り出したいというのは、タイル的な切り出しを行いたいという意図であることも多いと思います。このような時に地理座標を使う（＝ポリゴンで抽出）してしまうと、今回の例の「35度以上の範囲を切り出したつもりなのに 35.001度が含まれない」といったことが起こります。この場合、単純に緯度経度の数字だけで絞り込むほうが、目的に合致したデータが取得できるかもしれません。
どこまで正確性が要求されるのか、切り出したい目的は何なのかなど、こんな簡単な「四角形の範囲を切り出す」という処理ひとつとっても、考えることがたくさんあるのですね。

追記：PostGISでは挙動が異なる

　PostGISでは、このエントリで書いたMySQLでの挙動とは異なる結果を返しました。

テーブルとデータ作成（緯度経度がMySQLとは逆であることに注意）

CREATE TABLE points (p GEOMETRY(POINT, 4326));

INSERT INTO points VALUES (ST_GeomFromText('POINT(135.5 35.005)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(135.5 35.004)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(135.5 35.003)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(135.5 35.002)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(135.5 35.001)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(135.5 35)',4326));
INSERT INTO points VALUES (ST_GeomFromText('POINT(135.5 34.9999)',4326));

数値範囲での検索。当然、MySQLと同じ結果。

test=#  SELECT ST_AsText(p) FROM points
  WHERE ST_Y(p) BETWEEN 35 AND 36
    AND ST_X(p) BETWEEN 135 AND 136;
      st_astext      
---------------------
 POINT(135.5 35.005)
 POINT(135.5 35.004)
 POINT(135.5 35.003)
 POINT(135.5 35.002)
 POINT(135.5 35.001)
 POINT(135.5 35)

Polygonに内包する点の判定。MySQLには含まれなかった 35.001 が含まれている。35ジャストが含まれないのは、おそらく線上の点はWITHINではないのだろう。

test=# SELECT ST_AsText(p) FROM points 
  WHERE ST_WITHIN(p, ST_GeomFromText('POLYGON((135 35, 135 36, 136 36 ,136  35 , 135  35))',4326));
      st_astext      
---------------------
 POINT(135.5 35.005)
 POINT(135.5 35.004)
 POINT(135.5 35.003)
 POINT(135.5 35.002)
 POINT(135.5 35.001)

　この挙動について井口さんからコメントをいただきました（実際には順序が逆で、このコメントをいただいたから上記追加検証をしたのですが）。