MySQL GIS拡張関数： STX_Angle()

STX_Angle(p1, p2, p3)

　3つのPOINTを与えてそれが為す角（p1-p2-p3)の角度を求めます。戻り値はラジアン。PostGISに倣い、本関数は与えられた座標を平面として扱います。また、PostGISに実装されている「4つのPOINTを与える機能」「2つのLINESTRINGを与える機能」は実装されていません（要望があったら考えます）。
　角度はp1からp2を経由してp3に向かう p2における角度を、「反時計回りで」の角度として返します。

動作紹介

• 3時の方向から来て真上に抜けていく3点：反時計回りに90度

mysql> SELECT STX_Angle( ST_GeomFromText('POINT(10 10)'), 
    ->                   ST_GeomFromText('POINT( 5 10)'),
    ->                   ST_GeomFromText('POINT( 5 15)')) / PI()*180 agl; 
+------+
| agl  |
+------+
|   90 |
+------+
1 row in set (0.000 sec)

• 上から来て8時のほうに抜けていく線（左回りに90度をちょっと超えたくらい）

mysql> SELECT STX_Angle( ST_GeomFromText('POINT(15 15)'), 
    ->                   ST_GeomFromText('POINT(15 10)'),
    ->                   ST_GeomFromText('POINT( 5 5)')) / PI()*180 agl; 
+--------------------+
| agl                |
+--------------------+
| 116.56505117707799 |
+--------------------+
1 row in set (0.000 sec)

• 1時のほうから来て真上に抜ける角は38度くらい。STX_Azimuthと異なり、地球の丸さを考慮していない点に注意

mysql> SELECT STX_Angle( ST_GeomFromText('POINT(36.103056 140.088611)', 4326),
    ->                   ST_GeomFromText('POINT(35.658099 139.741357)', 6668),
    ->                   ST_GeomFromText('POINT(37        139.741357)', 6668))/ PI()*180 agl;
+-------------------+
| agl               |
+-------------------+
| 37.96924208395343 |
+-------------------+
1 row in set (0.000 sec)

• 参考：地球の丸さを考慮したSTX_Azimuthの結果。32.3度くらいになる。距離が長い（数十キロ）こともあり随分とSTX_Angleの結果とは異なります

mysql> SELECT STX_AZimuth(ST_GeomFromText('POINT(35.658099 139.741357)', 6668),
    ->                    ST_GeomFromText('POINT(36.103056 140.088611)', 4326))/ PI()*180 agl;
+-------------------+
| agl               |
+-------------------+
| 32.31909836897371 |
+-------------------+
1 row in set (0.000 sec)

補足説明

　冒頭でも解説した通り、この関数は地球の丸さを考慮しません。これはPostGISの実装に倣ったものです。この関数自体はOGCなどで定義されているものではなく、PostGIS独自開発のため、明確な共通仕様がないことから、MySQLでも同じ動作となるようにしました。とは言え、引数については今のところ 3つのPOINTを与えるのみ（PostGISは2つのLINESTRINGなどを受け取る仕様もある）なので、追加実装の必要が生じたら実装しようかなと思います。
　STX_Azimuthが北方向からターゲットPOINT方向への「時計回りの」角度を返すのに対して、本関数は p1からp3に向けての「反時計回りの」角度を返すことにも注意が必要です。

　また、STX_Azimuthとの差を見ると分かるとおり、地球規模の地点間の角度を見るのには向いていません。町の中など、ごく近い範囲の角度を得るために（平面での近似のため）高速な処理が得られるものとして割り切って使う関数かなと思います。

（このへんの仕様の違いを把握したり、それで良いのか作り替えるべきなのかを検討したりで、地味な関数ながら、ここまで書いてきた関数紹介の中で一番時間がかかった・・・・(日曜の午前が溶けた...)）

　余談ですが、地理座標系での動作テストのSQLと結果をAIさんに渡したら、「この東京から筑波への方位角は」と しれっと言っていて、たまげました。なんかもう、「察する能力が高い」と評価せざるを得ません。そう、日本経緯度原点からつくばVLBI(跡)への原点方位角なのです。