経度一度はどれくらいの長さ？MySQLをつかって調べてみよう

この日記は、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くらい）を歩いて体験できるしかけがあります。なかなか面白いですよ。

f:id:sakaik:20181220115619j:plain