🌀
Goによる台風経路図作成ロジック | Resilire Tech Blog

2024/11/13に公開
ポリゴン
 はじめにレジリア（Resilire） のバックエンドエンジニアKeita Miyanoです。
レジリアではサプライチェーンのリスク管理を行えるSaaSを提供しており、自然災害の影響を受けている可能性があるサプライチェーン上の拠点を可視化する機能があります。
災害情報の可視化ステップについてはこちらの記事も併せてご確認ください。
https://zenn.dev/resilire/articles/202410-using-jma-data

 台風経路図リスク可視化のひとつとして台風経路描画の機能を開発しています。毎年7-10月あたりは台風の発生が多く、この時期はニュースなどでこの経路図をよく見かけることがあると思います。（台風経路図の見方は、画像右下の通り）
この経路図は、気象庁が3時間ごとに発表している「台風解析・予報情報（延長予報）電文（新形式）」というxmlデータを元に描画できます。（こちらのURLから、VPTW60~VPTW65のリンクにアクセスすることで取得できます）
ただ、基本的に「n時間後における台風の中心と半径」しか描画に関係する数値は載ってません。（以下xmlの一部）
...

<ProbabilityCircle type="予報円">
    <jmx_eb:BasePoint description="北緯２６．８度東経１３７．８度" type="中心位置（度）"> +26.8+137.8/ </jmx_eb:BasePoint>
    <jmx_eb:BasePoint description="北緯２６度５０分東経１３７度５０分" type="中心位置（度分）"> +2650+13750/ </jmx_eb:BasePoint>
    <jmx_eb:Axes>
        <jmx_eb:Axis>
            <jmx_eb:Direction condition="全域" description="全域" type="方向" unit="８方位漢字"/>
            <jmx_eb:Radius description="中心が７０パーセントの確率で入る予報円の半径５７海里" type="７０パーセント確率半径" unit="海里"> 57 </jmx_eb:Radius>
            <jmx_eb:Radius description="中心が７０パーセントの確率で入る予報円の半径１０５キロ" type="７０パーセント確率半径" unit="km"> 105 </jmx_eb:Radius>
        </jmx_eb:Axis>
    </jmx_eb:Axes>
</ProbabilityCircle>

...
これだけの情報で暴風域を示す赤枠線などつくれる気がしなかったのですが、調べてもほぼ該当の記事が出てきませんでした。

そのため気象庁が公開しているロジック（※後述）を参考に、Goで台風経路図をポリゴン化するロジックを作成しました。
この記事では、GeoJSONファイルとして台風経路図を書き出すところまでをゴールとします。

 気象庁の描画ロジック気象庁による台風経路図の描画ロジックは公開されていて、配信資料に関する技術情報（気象編）第237号 というPDFファイルに詳細のロジックとFortranのサンプルコードが載っています。
左の暴風警戒域の円たちから、右の枠線を算出するロジックです。
見ていただければわかる通り、気象庁は非常に厳密な計算をしており場合分けの数も非常に多いです。（受験数学を思い出す）
気象庁の算出方針を平易にまとめると、
円同士を接線で結ぶ
内側にある不要な線を削除する
になるかと思います。
1.は、接線を引けない場合（片方の円がもう片方に含まれるなど）についての考慮をしたり、2.についても「不要な線はどこか？」ということに場合分けが必要になってきます。
しかし、最終的にポリゴン（点の集合）さえ得られればいいのでもう少し処理を簡素にできそうです。

 簡素なロジック暴風域の赤枠線のみつくることができれば、予報円も同様に作成できるため、ここでは暴風域に限定して書いていきます。
大きく分けて次の2つの手法でポリゴンを作っていきます。
凸包（convex hull）→ 円同士を接線で結ぶ
ポリゴンの自己交差削除 → 内側にある不要な線を削除する
簡単な数学と、地理空間の幾何ライブラリの機能を使って楽しようという考えです。

 凸包（convex hull）凸包とは与えられた集合を含む最小の凸集合であり、2次元においては与えられた点をすべて包含する最小の凸多角形のことです。
こちらのアルゴリズムをGoでConvexHull関数として実装すると以下のようになります。
import (
	"math"
	"sort"
)

type Point struct {
	Latitude  float64
	Longitude float64
}

// Helper function to determine the orientation of three points
// 0 -> p, q and r are collinear
// 1 -> Clockwise
// -1 -> Counterclockwise
func orientation(p, q, r Point) int {
	val := (q.Longitude-p.Longitude)*(r.Latitude-q.Latitude) - (q.Latitude-p.Latitude)*(r.Longitude-q.Longitude)
	if val == 0 {
		return 0
	}
	if val > 0 {
		return 1
	}
	return -1
}

// Distance between two points (Euclidean distance)
func dist(p1, p2 Point) float64 {
	return math.Sqrt((p2.Latitude-p1.Latitude)*(p2.Latitude-p1.Latitude) + (p2.Longitude-p1.Longitude)*(p2.Longitude-p1.Longitude))
}

func ConvexHull(points []Point) []Point {
	n := len(points)
	if n < 3 {
		return points
	}

	// Step 1: Find the point with the lowest Latitude (in case of tie, the leftmost Longitude)
	p0 := points[0]
	for i := 1; i < n; i++ {
		if points[i].Latitude < p0.Latitude || (points[i].Latitude == p0.Latitude && points[i].Longitude < p0.Longitude) {
			p0 = points[i]
		}
	}

	// Step 2: Sort the points based on polar angle with p0
	sort.Slice(points, func(i, j int) bool {
		o := orientation(p0, points[i], points[j])
		if o == 0 {
			return dist(p0, points[i]) < dist(p0, points[j])
		}
		return o == -1
	})

	// Step 3: Build the convex hull using a stack
	hull := []Point{p0, points[1], points[2]}

	for i := 3; i < n; i++ {
		for len(hull) > 1 && orientation(hull[len(hull)-2], hull[len(hull)-1], points[i]) != -1 {
			hull = hull[:len(hull)-1]
		}
		hull = append(hull, points[i])
	}

	return hull
}
凸包を利用して、2つの円を接線で結ぶことができます。
2つの円を構成するすべての点を対象として凸包を計算すると、下の画像のように内部の点群を除去して新たな1つの点群になります。どちらか一方が片方に包含されている場合には大きい方の円のみ残るため、場合分けのいらない手法です。
2つの円の点群
2つの円の点群へ凸包を適用
この性質を用いて、隣合う円のペアに対して凸包を計算すると以下のようになります。
次に、暴風域の枠内の線を除去していきます。

 ポリゴンの自己交差削除厳密にはポリゴンの自己交差修正です。

 GEOSのインストールこちらの処理にはGEOSというC++ベースのライブラリを使用するので、事前にインストールが必要です。
MacOSでは以下でインストールできます。
brew install geos
またGEOSのインターフェースであるgo-geosもインストールします。
go get github.com/twpayne/go-geos

 Bufferメソッド自己交差削除にBufferメソッドというものを使っていきます。
先ほど隣合うペアを凸包した後のポリゴンたちですが、一旦これらをマルチポリゴン（ポリゴンの集合体）として扱います。（後述のMultiPolygonToWKT関数で、以下のようなWKTという形式にします）
MULTIPOLYGON(
    ((139.70 35.65, 139.75 35.65, 139.75 35.70, 139.70 35.70, 139.70 35.65)),
    ((139.76 35.66, 139.81 35.66, 139.81 35.71, 139.76 35.71, 139.76 35.66))
)
ここでマルチポリゴンに対してGEOSのBufferメソッドを適用するのですが、本来こちらのメソッドは指定した距離でジオメトリの周囲にバッファを作成するためのメソッドです。（例えば、太さのない線に太さを付与できる）
// Buffer returns g with the given buffer.
func (g *Geom) Buffer(width float64, quadsegs int) *Geom {
	g.mustNotBeDestroyed()
	g.context.Lock()
	defer g.context.Unlock()
	return g.context.newNonNilGeom(C.GEOSBuffer_r(g.context.handle, g.geom, C.double(width), C.int(quadsegs)), nil)
}
第一引数のwidthはバッファの距離、第二引数のquadsegsは曲線をどれくらい滑らかにするか制御するパラメータです。（値が大きいほど滑らかになる）
このBufferメソッドですが、第一引数のwidthの値を0にするとポリゴンの自己交差を自動で修正し、今回の場合だと枠内の線を削除してくれます。
凸包で作られた円のペアに対して、これらの処理を実装すると以下のようになります。
import (
	"fmt"
	"strconv"
	"strings"

	"github.com/twpayne/go-geos"
)

func PointsToPolygonWKT(points []Point) string {
	// ポリゴンが閉じているかチェックし、閉じていない場合は閉じる
	if len(points) > 1 && (points[0].Latitude != points[len(points)-1].Latitude || points[0].Longitude != points[len(points)-1].Longitude) {
		points = append(points, points[0])
	}

	var coords []string
	for _, point := range points {
		coords = append(coords, fmt.Sprintf("%f %f", point.Longitude, point.Latitude))
	}

	// POLYGON WKT形式に変換
	wkt := fmt.Sprintf("((%s))", strings.Join(coords, ", "))
	return wkt
}

// [][]PointからWKT形式のMULTIPOLYGONを作成する関数
func MultiPolygonToWKT(multiPolygon [][]Point) string {
	var polygons []string
	for _, polygon := range multiPolygon {
		polygons = append(polygons, PointsToPolygonWKT(polygon))
	}

	// MULTIPOLYGON WKT形式に変換
	wkt := fmt.Sprintf("MULTIPOLYGON(%s)", strings.Join(polygons, ", "))
	return wkt
}

func WKTToPolygonPoints(wkt string) ([]Point, error) {
	// WKTからPOLYGONの座標部分を抽出
	wkt = strings.TrimPrefix(wkt, "POLYGON ((")
	wkt = strings.TrimSuffix(wkt, "))")
	// まれに内側にポリゴンが残ることがあるので削除
	wkt = strings.Split(wkt, "), (")[0]
	coordPairs := strings.Split(wkt, ", ")

	var polygon []Point

	// 各座標ペアを処理
	for _, pair := range coordPairs {
		coords := strings.Split(pair, " ")
		if len(coords) != 2 {
			// 詳細なエラーメッセージを出力
			fmt.Printf("無効な座標ペア: %s\n", pair)
			continue // 無効な座標ペアをスキップ
		}

		// 経度と緯度をfloat64に変換
		lon, err := strconv.ParseFloat(coords[0], 64)
		if err != nil {
			return nil, fmt.Errorf("無効な経度値: %v", err)
		}
		lat, err := strconv.ParseFloat(coords[1], 64)
		if err != nil {
			return nil, fmt.Errorf("無効な緯度値: %v", err)
		}

		// Point構造体に追加
		polygon = append(polygon, Point{Latitude: lat, Longitude: lon})
	}

	return polygon, nil
}

func RemoveSelfIntersections(stormAreaPairs [][]Point) []Point {
	wkt := MultiPolygonToWKT(stormAreaPairs)
	geom, err := geos.NewGeomFromWKT(wkt)
	if err != nil {
		log.Fatalf("CalcStormAreaPolygon Error: %v, wkt: %v", err, wkt)
	}
	buffered := geom.Buffer(0, 32)

	bufferedWKT := buffered.ToWKT()

	bufferedPoints, err := WKTToPolygonPoints(bufferedWKT)
	if err != nil {
		log.Fatalf("CalcStormAreaPolygon Error: %v, polygon: %v", err, bufferedWKT)
	}

	return bufferedPoints
}
実際にここまでの処理を実行すると、枠内の線が消えて枠線のみが残ります。

 GeoJSONファイル書き出しGeoJSONの書き出しにはgo.geojsonを使用するので、こちらもインストールが必要です。
go get github.com/paulmach/go.geojson
GeoJSON書き出しは本題ではないので説明は省きますが、点群ができてしまえば書き出しは簡単に行えます。今回はgeojson.ioで手軽に確認することを目的としているため、SetPropertyで線の色や太さを指定しています。
import (
	"encoding/json"
	"fmt"
	"os"

	geojson "github.com/paulmach/go.geojson"
)

func MakeGeojsonPolygon(points []Point) *geojson.Feature {
	geoJsonPoints := make([][]float64, 0, len(points)+1)
	for _, coordinate := range points {
		geoJsonPoints = append(
			geoJsonPoints,
			[]float64{coordinate.Longitude, coordinate.Latitude},
		)
	}
	// 最初の点と最後の点を一致させる
	geoJsonPoints = append(
		geoJsonPoints,
		[]float64{points[0].Longitude, points[0].Latitude},
	)
	coordinates := [][][]float64{geoJsonPoints}
	polygon := geojson.NewPolygonFeature(coordinates)
	return polygon
}

func SaveGeoJSONToFile(filename string, data []byte) error {
	return os.WriteFile(filename, data, 0644)
}

func WriteGeoJSON(stormAreaBorderPoints []Point) {
	featureCollection := geojson.NewFeatureCollection()
	stormAreaBorderPolygon := MakeGeojsonPolygon(stormAreaBorderPoints)
	// 線の色や太さ指定
	stormAreaBorderPolygon.SetProperty("stroke", "#ff0000")
	stormAreaBorderPolygon.SetProperty("stroke-width", 2)
	stormAreaBorderPolygon.SetProperty("stroke-opacity", 1)
	stormAreaBorderPolygon.SetProperty("fill-opacity", 0)
	featureCollection.AddFeature(stormAreaBorderPolygon)

	// GeoJSONとしてエンコード
	geoJSON, err := json.MarshalIndent(featureCollection, "", "  ")
	if err != nil {
		fmt.Println("Error encoding GeoJSON:", err)
		return
	}

	// ファイルに保存する
	savePath := "output.geojson"
	err = SaveGeoJSONToFile(savePath, geoJSON)
	if err != nil {
		fmt.Println("Error saving GeoJSON to file:", err)
		return
	}

	fmt.Printf("GeoJSON successfully written to %s\n", savePath)
}
算出したポリゴンをGeoJSON化し、geojson.ioで表示させた結果が以下です。（暴風域と予報円、中心線のみ）
2024年8月9月の台風データに対して今回のロジックを適用してみましたが、問題なくポリゴンが算出されました。
ちなみに1つの台風経路図に対してGeoJSONファイルの書き出しにかかる時間は、私の環境で10~20msほどでしたので実用的です。事前に計算したポリゴンをどこかに格納することなく、リアルタイムに計算できるかと思います。

 おわりにGoによる台風経路図の作成について書きました。台風経路図そのものはさまざまなサイトで見ることはできるものの、サプライチェーン上の拠点とともに表示することで影響予測などを効率的に行うことができます。
気象庁データだけでなく、災害情報は他にも地図に可視化することで多くの示唆を生むデータが多く、奥の深い領域です。

 エンジニア採用強化中レジリアでは、サプライチェーンリスク管理クラウドサービスResilireの開発メンバーを積極募集中です！
https://recruit.resilire.jp/