気象庁XMLの台風解析・予報情報電文を使用して台風情報を描画する備忘録 その2

この記事は、防災アプリ Advent Calendar 2023 20日目の記事です。
[気象庁XMLの台風解析・予報情報電文を使用して台風情報を描画する備忘録 その1](https://qiita.com/oruponu/items/2c1637c226bc9a61b43f) の続きです。


# 前回のあらすじ
前回は台風の現在位置（強風域、暴風域）を表示するところまでを記事にしました。
今回は台風の進路予報を描画するところまでを書きます。


# 予報円を描画

![46d6d1fb-db7c-4dd6-a956-ba05574f7750.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/195741/bdaed535-79d7-566a-f48c-87874ad0f49e.png)

`Body/MeteorologicalInfos[0]/MeteorologicalInfo` が配列になっており、この中に予報の時間ごとの要素が格納されています。
この配列の各要素に対してループ処理を行います。

まず、`Body/Meteorologic
alInfos[0]/MeteorologicalInfo[i]/Item[0]/Kind` から `Property/Type` が「中心」の要素を取得します。
この要素の中から `Property[0]` を変数 `center` に置いておきます。

---

次に、`center` の中に予報円の要素があるか確認します。
予報円の要素は `center/CenterPart/ProbabilityCircle[@type='予報円']` で表されます。

予報円の要素があれば、強風域や暴風域を描画したときのように `center/CenterPart/ProbabilityCircle[0]/BasePoint[@type='中心位置（度）']` を中心とし、半径が `center/CenterPart/ProbabilityCircle[0]/Axes/Axis[0]/Radius[@unit='km']` km の円を描画します。


# 共通外接線を描画

![1e2c37d3-a43c-4c7b-9c93-69332934e5e4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/195741/ca2747bc-cc81-1e96-06f5-c3f52b840a2b.png)

共通外接線の線分を描画するには、線分の始点の円と終点の円のそれぞれの位置を用いる必要があります。ループ処理をするために、予報円の位置と半径をあらかじめ配列に保存しておきます。
さらに、台風の中心から1つ目の予報円に向かう線も描画する必要があります。これは、台風の中心位置に半径0の円があると仮定して、この情報も配列に追加しておきます。

今回は共通外接線の計算を簡単にするため、**円の中心座標を緯度経度からピクセル座標に変換してから配列に追加しました**。[^1]
円の中心座標（ピクセル座標）と円の半径を保存した配列を `forecastCircles` としておきます。

---

`forecastCircles` に対して `forecastCircles` の要素数 - 1回だけループを回します。ループの中で隣り合う2つの円の情報を取得します。
2つの円に対して、共通外接線を計算する関数を呼び出します。

以下は、共通外接線を計算する JavaScript のサンプルコードです。
```javascript
/*
 * x1: 1つ目の円の中心のX座標
 * y1: 1つ目の円の中心のY座標
 * r1: 1つ目の円の半径
 * x2: 2つ目の円の中心のX座標
 * y2: 2つ目の円の中心のY座標
 * r2: 2つ目の円の半径
 */
function calcOuterTangent(x1, y1, r1, x2, y2, r2) {
  const dx = x2 - x1;
  const dy = y2 - y1;
  const dist = Math.sqrt(dx * dx + dy * dy);

  // 一方の円がもう一方の円の内部にある場合は何も返さない
  if (dist <= Math.abs(r2 - r1)) return;

  const angle1 = Math.atan2(dy, dx);
  const angle2 = Math.acos((r1 - r2) / dist);

  const calcTangentPoint = (x, y, radius, offset) => [
    x + radius * Math.cos(angle1 + offset),
    y + radius * Math.sin(angle1 + offset)
  ];

  return [
    [
      calcTangentPoint(x1, y1, r1, angle2),
      calcTangentPoint(x2, y2, r2, angle2)
    ],
    [
      calcTangentPoint(x1, y1, r1, -angle2),
      calcTangentPoint(x2, y2, r2, -angle2)
    ]
  ];
}
```

戻り値は `outerTangentPoints` としておきます。
一方の円がもう一方の円の内部にある場合、共通外接線は存在しないため描画しません。その状態のとき `outerTangentPoints` は `undefined` となっています。

`outerTangentPoints` は共通外接線の端点の座標を表す配列で、次のような構造になっています。

- 1本目の線分：
    - X1：`outerTangentPoints[0][0][0]`
    - Y1：`outerTangentPoints[0][0][1]`
    - X2：`outerTangentPoints[0][1][0]`
    - Y2：`outerTangentPoints[0][1][1]`
- 2本目の線分：
    - X1：`outerTangentPoints[1][0][0]`
    - Y1：`outerTangentPoints[1][0][1]`
    - X2：`outerTangentPoints[1][1][0]`
    - Y2：`outerTangentPoints[1][1][1]`


# 暴風警戒域を描画（概要のみ）

![ad8b39c3-bd41-4841-a4dd-7ec8286a8fef.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/195741/ad1697a2-2779-2aa1-8b6d-375b7c9742e9.png)

詳しく説明すると長くなりそうだったので、本記事では処理の概要だけ簡単に記載します。

処理は大きく3ステップに分かれます。
|1. 円の座標を計算|2. 共通外接線の座標を計算|3. ポリゴンをマージ|
|---|---|---|
|![e47667d3-9829-f473-d303-d31b6c17b266.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/195741/e47667d3-9829-f473-d303-d31b6c17b266.png)|![680b910e-7323-1476-5d4b-d6ea4569da38.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/195741/680b910e-7323-1476-5d4b-d6ea4569da38.png)|![a81a80e9-f142-72c2-c63c-c86f4337c834.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/195741/a81a80e9-f142-72c2-c63c-c86f4337c834.png)|

1と2については、予報円とほとんど同様の処理を行います。
予報円と違うのは、暴風域も一つの円として扱う必要があることと、暴風警戒域がある時間帯の間に暴風警戒域のない時間帯が挟まる可能性を考慮する必要があることです。

共通外接線は線分同士を結んで閉じた図形にしておきます。

<img width=250 src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/195741/6e2652c7-d39e-b97d-0465-5e97f80fec04.png>


最後に、円と共通外接線を含む四角形をマージして[^2]描画します。

# おわりに
本記事では、台風の進路予報を描画するところまでを紹介しました。
もっと効率的な手法はあるかもしれませんが、一つのアプローチとして参考程度にご覧いただければ幸いです。

また、実際の電文を受け取って描画する場合は、さらに考慮しないといけないこともあります。
予報円については、通常24時間刻みの予報が発表されますが、台風が日本に接近して台風情報の発表間隔が1時間ごとになると24時間先までの予報が3時間刻みになります。
そのまま描画すると予報円が密集して視認性が損なわれることもあるため、必要に応じて間引く処理が必要です。

今回までの内容で取りこぼした内容もあるので、機会があれば続きを書こうと思います。

[^1]: サンプル画像のコードでは D3.js の投影関数を使用して変換しています
[^2]: サンプル画像のコードでは Turf.js の機能でマージしています