この日記は、RDBMS-GIS(MySQL,PostgreSQLなど) Advent Calendar 2019 の2日目の記事です。
( https://qiita.com/advent-calendar/2019/rdbms_gis )
経度一度ってどれくらい？
赤道付近の一周の長さが だいたい 40,000km だというのはみんな知っていますよね。実際はもう少し長いのですが、とりあえず感覚的なものでいいです。で、経度というのは、緯度が上がっていくごとに円が小さくなりますから、1度の長さもどんどん短くなるはずです。どんな感じなのかな、とMySQLを使って見てみました。
MySQL、まだまだ関数が少ないので（とくに集約系や分析的なもの）、関数が少ないうちに、少ないこれらの関数を遊び倒しましょう。できることが限られているというのは、アイデアを膨らませるという点でメリットでもあります！
経度1度を求める
まず、テーブル g1 を作って適当にデータを入れます。ここでは、緯度0度から10度刻みで、緯度135度と136度の2つの位置情報をカラムに入れることにしました。測地系は、JGD2011を使います。
CREATE TABLE g1 (
id integer,
pos1 GEOMETRY SRID 6668,
pos2 GEOMETRY SRID 6668
);
INSERT INTO g1 VALUES (0, ST_GeomFromText('POINT(0 135)', 6668), ST_GeomFromText('POINT(0 136)', 6668));
INSERT INTO g1 VALUES (10, ST_GeomFromText('POINT(10 135)', 6668), ST_GeomFromText('POINT(10 136)', 6668));
INSERT INTO g1 VALUES (20, ST_GeomFromText('POINT(20 135)', 6668), ST_GeomFromText('POINT(20 136)', 6668));
INSERT INTO g1 VALUES (30, ST_GeomFromText('POINT(30 135)', 6668), ST_GeomFromText('POINT(30 136)', 6668));
INSERT INTO g1 VALUES (35, ST_GeomFromText('POINT(35 135)', 6668), ST_GeomFromText('POINT(35 136)', 6668));
INSERT INTO g1 VALUES (40, ST_GeomFromText('POINT(40 135)', 6668), ST_GeomFromText('POINT(40 136)', 6668));
INSERT INTO g1 VALUES (50, ST_GeomFromText('POINT(50 135)', 6668), ST_GeomFromText('POINT(50 136)', 6668));
INSERT INTO g1 VALUES (60, ST_GeomFromText('POINT(60 135)', 6668), ST_GeomFromText('POINT(60 136)', 6668));
INSERT INTO g1 VALUES (70, ST_GeomFromText('POINT(70 135)', 6668), ST_GeomFromText('POINT(70 136)', 6668));
INSERT INTO g1 VALUES (80, ST_GeomFromText('POINT(80 135)', 6668), ST_GeomFromText('POINT(80 136)', 6668));
INSERT INTO g1 VALUES (90, ST_GeomFromText('POINT(90 135)', 6668), ST_GeomFromText('POINT(90 136)', 6668));
このテーブルの各レコードの長さは、以下のSQLで求めることができます。
SELECT id, ST_Distance(pos1, pos2) FROM g1 ORDER BY ID;
+------+-------------------------+
| id | ST_Distance(pos1, pos2) |
+------+-------------------------+
| 0 | 111319.49079326246 |
| 10 | 109639.3390102644 |
| 20 | 104646.97563995633 |
| 30 | 96486.0081494232 |
| 35 | 91287.79089978167 |
| 40 | 85393.3619519223 |
| 50 | 71695.04015879167 |
| 60 | 55799.17666195067 |
| 70 | 38185.80088710373 |
| 80 | 19393.044953483724|
| 90 | 0 |
+------+-------------------------+
11 rows in set (0.00 sec)
赤道付近で、1度あたり111km、日本がある北緯35度付近では91km、70度で38km、80度で19kmと一気に狭くなっていくのがわかりますね。感覚とも一致します。
1度ではなく1分は？
何十キロもの長さは日常でイメージしにくいので、もう少し小さい単位である「分」の単位で距離を求めてみましょう。1分は約 0.0167度 なので以下のようなデータを作ります。各緯度で、経度135度から135.0167度を表すデータ群です。
CREATE TABLE g2 (id integer, pos1 GEOMETRY SRID 6668, pos2 GEOMETRY SRID 6668);
INSERT INTO g2 VALUES (0, ST_GeomFromText('POINT(0 135)', 6668), ST_GeomFromText('POINT(0 135.0167)', 6668));
INSERT INTO g2 VALUES (10, ST_GeomFromText('POINT(10 135)', 6668), ST_GeomFromText('POINT(10 135.0167)', 6668));
INSERT INTO g2 VALUES (20, ST_GeomFromText('POINT(20 135)', 6668), ST_GeomFromText('POINT(20 135.0167)', 6668));
INSERT INTO g2 VALUES (30, ST_GeomFromText('POINT(30 135)', 6668), ST_GeomFromText('POINT(30 135.0167)', 6668));
INSERT INTO g2 VALUES (35, ST_GeomFromText('POINT(35 135)', 6668), ST_GeomFromText('POINT(35 135.0167)', 6668));
INSERT INTO g2 VALUES (40, ST_GeomFromText('POINT(40 135)', 6668), ST_GeomFromText('POINT(40 135.0167)', 6668));
INSERT INTO g2 VALUES (50, ST_GeomFromText('POINT(50 135)', 6668), ST_GeomFromText('POINT(50 135.0167)', 6668));
INSERT INTO g2 VALUES (60, ST_GeomFromText('POINT(60 135)', 6668), ST_GeomFromText('POINT(60 135.0167)', 6668));
INSERT INTO g2 VALUES (70, ST_GeomFromText('POINT(70 135)', 6668), ST_GeomFromText('POINT(70 135.0167)', 6668));
INSERT INTO g2 VALUES (80, ST_GeomFromText('POINT(80 135)', 6668), ST_GeomFromText('POINT(80 135.0167)', 6668));
INSERT INTO g2 VALUES (90, ST_GeomFromText('POINT(90 135)', 6668), ST_GeomFromText('POINT(90 135.0167)', 6668));
SELECT id, ST_Distance(pos1, pos2) FROM g2 ORDER BY ID;
+------+-------------------------+
| id | ST_Distance(pos1, pos2) |
+------+-------------------------+
| 0 | 1859.0354970687047 |
| 10 | 1830.9776594574585 |
| 20 | 1747.607087970476 |
| 30 | 1611.3214483744373 |
| 35 | 1524.512474129115 |
| 40 | 1426.0766207422205 |
| 50 | 1197.3160879967666 |
| 60 | 931.8551188371557 |
| 70 | 637.7100204965047 |
| 80 | 323.86784178140675|
| 90 | 0 |
+------+-------------------------+
赤道付近で1.9km、北緯35度付近で1.5kmです。念のため赤道付近の「分」の結果を60倍して答え合わせをしてみましょうか。
mysql> SELECT 1859.0354970687047 * 60 ;
+-----------------------+
| 1859.0354970687047*60 |
+-----------------------+
| 111542.1298241222820 |
+-----------------------+
1 row in set (0.00 sec)
あれれ。本来の赤道付近の距離 111319.49079326246 よりも 少し大きくなっていまいましたね。考えてみたら、
「1度」を 0.0167 としましたが、これがちょっと切り上げが雑すぎたのかもしれません。0.016666666666667くらいでやるべきだったのでしょうか。まぁこれは感覚とも一致するのでよしとしましょう。
というかそもそも、0.0166666....とか言っている時点で誤差許容なので、もうちょっとまともな計算をしましょうか。
経度1度、1分、1秒の長さを求める(別解)
最初に「度」を求めたのですから、度を60で割って分を求めたり、度を60*60で割って秒を求めたりするほうが(先ほどみたいに途中で不正確な小数を持ち出すよりも）正確な値が求まります。ここでは、そのやり方で。最初に作った g1 テーブルを使います。
WITH cte_g1 AS (
SELECT id, ST_Distance(pos1, pos2) dist FROM g1
)
SELECT id, round(dist,2) do, round(dist/60,2) fun, round(dist/60/60,2) byou
FROM cte_g1;
+------+-----------+---------+-------+
| id | do | fun | byou |
+------+-----------+---------+-------+
| 0 | 111319.49 | 1855.32 | 30.92 |
| 10 | 109639.34 | 1827.32 | 30.46 |
| 20 | 104646.98 | 1744.12 | 29.07 |
| 30 | 96486.01 | 1608.10 | 26.80 |
| 35 | 91287.79 | 1521.46 | 25.36 |
| 40 | 85393.36 | 1423.22 | 23.72 |
| 50 | 71695.04 | 1194.92 | 19.92 |
| 60 | 55799.18 | 929.99 | 15.50 |
| 70 | 38185.80 | 636.43 | 10.61 |
| 80 | 19393.04 | 323.22 | 5.39 |
| 90 | 0.00 | 0.00 | 0.00 |
+------+-----------+---------+-------+
11 rows in set (0.00 sec)
うーん、CTE便利！ MySQL 8.0 最高！ まぁこの程度なら、FROMのサブクエリに書いてもいいんですけどね。
少数以下が長くあっても仕方がないので、ここでは小数第2位（センチメートル）で丸めて出力しました。
赤道付近での1秒は30.92m、北緯35度では 25.36m だとわかります。これくらいの長さだと、そのへんの通り道を眺めながら「あそこまでで1秒かー」とイメージが湧きやすいですね。
緯度はどうなるか
ここまでやったらならば、一方の緯度はどうなるか確認してみたくなります。
緯度は経度みたいには細くなっていかないので、一定のままだろう、と予想して、テーブルg3とそのデータを作って試します。
CREATE TABLE g3 (id integer, pos1 GEOMETRY SRID 6668, pos2 GEOMETRY SRID 6668);
INSERT INTO g3 VALUES (0, ST_GeomFromText('POINT(0 135)', 6668), ST_GeomFromText('POINT(1 135)', 6668));
INSERT INTO g3 VALUES (10, ST_GeomFromText('POINT(10 135)', 6668), ST_GeomFromText('POINT(11 135)', 6668));
INSERT INTO g3 VALUES (20, ST_GeomFromText('POINT(20 135)', 6668), ST_GeomFromText('POINT(21 135)', 6668));
INSERT INTO g3 VALUES (30, ST_GeomFromText('POINT(30 135)', 6668), ST_GeomFromText('POINT(31 135)', 6668));
INSERT INTO g3 VALUES (35, ST_GeomFromText('POINT(35 135)', 6668), ST_GeomFromText('POINT(36 135)', 6668));
INSERT INTO g3 VALUES (40, ST_GeomFromText('POINT(40 135)', 6668), ST_GeomFromText('POINT(41 135)', 6668));
INSERT INTO g3 VALUES (50, ST_GeomFromText('POINT(50 135)', 6668), ST_GeomFromText('POINT(51 135)', 6668));
INSERT INTO g3 VALUES (60, ST_GeomFromText('POINT(60 135)', 6668), ST_GeomFromText('POINT(61 135)', 6668));
INSERT INTO g3 VALUES (70, ST_GeomFromText('POINT(70 135)', 6668), ST_GeomFromText('POINT(71 135)', 6668));
INSERT INTO g3 VALUES (80, ST_GeomFromText('POINT(80 135)', 6668), ST_GeomFromText('POINT(81 135)', 6668));
INSERT INTO g3 VALUES (90, ST_GeomFromText('POINT(89 135)', 6668), ST_GeomFromText('POINT(90 135)', 6668));
mysql> SELECT id, ST_Distance(pos1, pos2) FROM g3 ORDER BY ID;
+------+-------------------------+
| id | ST_Distance(pos1, pos2) |
+------+-------------------------+
| 0 | 110573.13812416371 |
| 10 | 110610.23595555034 |
| 20 | 110710.37167306663 |
| 30 | 110861.46743133431 |
| 35 | 110950.61420268985 |
| 40 | 111045.29885197058 |
| 50 | 111239.69315258414 |
| 60 | 111421.20351156592 |
| 70 | 111567.93710081652 |
| 80 | 111662.1956839133 |
| 90 | 111692.61028462648 |
+------+-------------------------+
あれれ。みんなだいたい111kmではありますが（赤道付近の経度1度ともほぼ合致していますね）、よく見ると少しずつ値が異なりますね。そうか、地球（JGD2011でのモデルであるGRS80）は真球ではなく回転楕円体だから、緯度があがるごとに少しずつ値が変化するわけですね。確かにそうだ。
いや、でもよく見てみると、、、、緯度があがるほうが1度あたりの距離が長くなっていますね。イメージを真球から極半径を共有した回転楕円体へと拡張していくと、赤道付近のほうがぐーっと伸びているのだから、そちらのほうが1度あたりの距離が(真球と比べて)より長くなるようなイメージを持っていたのですが・・・・・かるく調べてみると値が間違っているわけではないようですが、感覚的に納得できません。。ここまでで今日の私はギブアップ。
誰か、「なぜ極付近のほうが緯度1度の距離が長くなるのか」、わかりやすくおしえてください！
おまけ
つくば市の国土地理院内にある「地図と測量の科学館」には、地球のミニチュアモデルがあります。曲面に乗っかった日本の形とか、なかなか面白いものですよ。そしてその周りには、「1秒の距離をあるいてみよう」という看板と標石が設置されていて、実際にその場所での1秒（南北30.8m、東西25mくらい）を歩いて体験できるしかけがあります。なかなか面白いですよ。
