気象庁が発表する台風予報を描画する際、各時刻の暴風警戒域を単に重ねて描画すると、円が重なって読み取りづらくなる場合があります。このため、暴風警戒域を表示する際には円と接線に囲まれる範囲を包絡する線で描画することが認められています
この暴風警戒域を包絡する範囲の計算は配信資料に関する技術情報(気象編)第237号(気象庁)に掲載されていますが、複雑な場合分けもあり、1500行ほどのプログラムとなっており、生成AIにJavaScriptに書き換えるようお願いしてみましたが断られてしまいました。トホホ
今回は、暴風警戒域を包絡する範囲の計算を簡易的に、ピュアJavaScriptで実装するよう試みます
なぜ導出手法が複雑になるか
気象庁の技術情報の導出手法は、暴風警戒域を包絡する範囲を計算する際に、まず円とその接線を計算し、その後接線がなす台形や円の内側の部分を除去するアルゴリズムで計算しています
ただ、この「内側にある」を判断する方法が難しいです。たとえば2つの円が接する位置関係は、円が他方の円にすっぽり収まる場合、円と円が一部重なり合う場合、円と円が離れて存在する場合が考えられますが、それぞれ消去される範囲が異なります(前者から 小さいほうの円全体、重なり合う部分の弧、消去される線分なし)
台形に関しては2点で交差するか4点で交差するか等さらに面倒です
気象庁の解説資料はこれらを真摯に考慮した結果、A4 12ページに及ぶ資料となっています。まるで大学入試の数学の解答を読んでいるかのような緻密さです
簡便に計算する手法を考える
図形と図形の交点を求めようとすると円と円なのか、円と直線なのか、内包関係は……など場合分けが複雑になるため、まず包絡線を構成する可能性のある点を多数求め、内側にある点を消去する方針とします
これによって、図形と図形の内包問題ではなく図形と点の内包問題となり、点が円の内側にあるか(距離と半径の比較で容易に判定可能)と点が台形の内側にあるか(気象庁の手法で比較的容易に判定可能)の2パターンのみの場合分けで済むようになります
大まかな方針
引数は各時刻の暴風警戒域の中心と半径、返り値は暴風警戒域を包絡する座標群とします
となりあう時刻の暴風警戒域の共通外接線と接点を求める
(これと同時に共通外接線の接点4点がなす台形が求まります)
↓
共通外接線と円を多数の点に分ける
↓
各点について、各時刻の円や台形の内側にあるか判定し、内側の場合は消去する
(円は中心と点の距離、半径を比較することで内包しているか判定)
(台形は気象庁の手法を参考に、三角形2つに分けて内包しているかを判定)
↓
(以上の手順でぶつ切りになった線分のうち、必要に応じ始点終点が近い線分どうしを接続する)
↓
点の集まりを返す
計算用のコード
以下がメインとなる関数です。circlesに暴風警戒域の円群([緯度, 経度, 半径(km)], [緯度, 経度, 半径(km)]...)を、resolutionKmでどの程度の精度で点をばらまくかを与えます
calcEggplantCoordinates()
function calcEggplantCoordinates( circles, resolutionKm=5){
const circleTrapezoidPoints = calcCircleTrapezoidPoints( circles, resolutionKm); // 円と接線上の多数の点の座標を計算
// return circleTrapezoidPoints.coordinates; // ここで返すとすべての円と接線が書かれた状態
const eggplantCoordinates = removeInnerPoints( circles, circleTrapezoidPoints.trapezoids, circleTrapezoidPoints.coordinates, resolutionKm); // 円と台形内の点を消去
// return eggplantCoordinates; // ここで返すとぶつ切りになった線分が返る
return jointCoordinates( eggplantCoordinates, resolutionKm); // 始点終点が近い線分を結合
}
円と接線上の多数の点の座標を求めます
calcCircleTrapezoidPoints()
function calcCircleTrapezoidPoints( circles, resolutionKm=5){
let trapezoids = [], coordinates = [];
for( let i=0; i<circles.length; i++){
let lat1 = circles[i][0], lon1 = circles[i][1], r1 = circles[i][2];
coordinates.push( generateCircleCoordinates( lat1, lon1, r1, resolutionKm)); // 円上の座標群を求める
if( i!=circles.length-1){ // 最後の時刻でない場合、次の時刻の間の接線を求める
let lat2 = circles[i+1][0], lon2 = circles[i+1][1], r2 = circles[i+1][2];
const tangents = getExternalTangents(lat1, lon1, r1, lat2, lon2, r2); // 接点の座標を求める
for( let coordinate of complementExternalTangents( tangents, resolutionKm)){ // 接点どうしを内分する点群を求める
coordinates.push( coordinate);
}
trapezoids.push( [tangents[0][0], tangents[0][1], tangents[1][1], tangents[1][0]]); // 内分判定のために台形の座標も必要なため、別に記録しておく
}
}
return { coordinates, trapezoids};
}
円上の各点の緯度経度を求めます。解像度に応じて何点に分けるかを決め、各方向の点の座標を求めていきます
// 特定の緯度経度を中心とした円上の各点の緯度・経度群を求める関数(中心緯度, 中心経度, 半径, 解像度)
function generateCircleCoordinates(centerLat, centerLon, radiusKm, resolutionKm = 10) {
const coords = [];
const points = Math.ceil( 2 * radiusKm * Math.PI / resolutionKm); // 解像度に応じて何点に分けるかを決める
for (let i = 0; i < points; i++) {
const angle = (i / points) * 360;
coords.push(getCirclePoint(centerLat, centerLon, radiusKm, angle));
}
return coords;
}
// 特定の緯度経度から特定の距離・角度にある点の緯度経度を求める関数(緯度, 経度, 距離, 角度)
function getCirclePoint(centerLat, centerLon, radiusKm, angleDeg) {
const R = 6378.137; // 地球の半径 (km)
const radLat = (centerLat * Math.PI) / 180, radLon = (centerLon * Math.PI) / 180;
const radAngle = (angleDeg * Math.PI) / 180;
const dR = radiusKm / R; // 角度としての距離
// 緯度の計算
const lat2 = Math.asin(
Math.sin(radLat) * Math.cos(dR) +
Math.cos(radLat) * Math.sin(dR) * Math.cos(radAngle)
);
// 経度の計算
const lon2 =
radLon +
Math.atan2(
Math.sin(radAngle) * Math.sin(dR) * Math.cos(radLat),
Math.cos(dR) - Math.sin(radLat) * Math.sin(lat2)
);
return {
lat: (lat2 * 180) / Math.PI,
lon: (lon2 * 180) / Math.PI,
};
}
接線上の各点の緯度経度を求めます
// 緯度経度をローカル直交座標(m)へ変換
function latLonToXY(lat, lon, lat0, lon0) {
// 中略。記事末尾のサンプルコードを参照
return { x, y};
}
// ローカル直交座標(m)を緯度経度へ変換
function xyToLatLon(x, y, lat0, lon0) {
// 中略。記事末尾のサンプルコードを参照
return { lat, lon};
}
// 2円の共通外接線の接点を求める。返り値は [ [[lat,lon], [lat,lon]], ... ]
function getExternalTangents( lat1, lon1, r1, lat2, lon2, r2){
// 中略。記事末尾のサンプルコードを参照
return result;
}
// 共通外接線の接点間を補完する関数
function complementExternalTangents( tangents0, resolutionKm=5) {
let tangents = [];
for( let tangent0 of tangents0){
const distanceKm = distance( tangent0[0].lat, tangent0[0].lon, tangent0[1].lat, tangent0[1].lon) / 1000;
const points = Math.ceil( distanceKm / resolutionKm); // 解像度に応じて何点に分けるかを決定
let tangent = [];
for (let i = 0; i < points; i++) {
tangent.push({ // 単純に始点終点を内分
lat: tangent0[0].lat + (tangent0[1].lat - tangent0[0].lat) * i / points,
lon: tangent0[0].lon + (tangent0[1].lon - tangent0[0].lon) * i / points
})
}
tangents.push( tangent);
}
return tangents;
}
// 緯度経度から距離を求める関数。返り値はメートル単位。https://zenn.dev/bellbind/scraps/bb2cfb915430f5 を参考に作成
function distance( lat1, lon1, lat2, lon2){
const rad = Math.PI / 180, r = 6378137;
const p1 = lat1 * rad, p2 = lat2 * rad, q1 = lon1 * rad, q2 = lon2 * rad;
return r * Math.acos(Math.cos(p1) * Math.cos(p2) * Math.cos(q1 - q2) + Math.sin(p1) * Math.sin(p2));
}
求めた座標群から円・台形の内部の点を除去します
function removeInnerPoints( circles, trapezoids, coordinates0, resolutionKm=5){
let coordinates = [], coordinate = [];
for( let coordinate0 of coordinates0){
for( let point of coordinate0){
let lat = point.lat, lon = point.lon, isInner = false;
for( let circle of circles){
if( distance( lat, lon, circle[0], circle[1])/1000 < ( circle[2] - resolutionKm/10) ){ // 点と円の内包単に距離で判定。解像度(resolutionKm)の10分の1だけ余裕をとる
isInner = true;
break;
}
}
for( let trapezoid of trapezoids){
if( isInnerTrapezoid( point, trapezoid, resolutionKm)){ // 点と台形の内包は専用の関数で判定
isInner = true;
break;
}
}
if( isInner){ // 内側にあると判断したら、それまでの座標群を一度返す
if( coordinate.length>0){
coordinates.push( coordinate);
}
coordinate = [];
}else{
coordinate.push( point);
}
}
if( coordinate.length>0){
coordinates.push( coordinate);
}
coordinate = [];
}
return coordinates;
}
点と台形の内包関係は、気象庁の技術情報掲載の通り、台形を三角形に分けて調べます
// 特定の点が台形内に存在するかを判定
function isInnerTrapezoid( point, trapezoid0, resolutionKm=5){
const gravity = { lat:(trapezoid0[0].lat+trapezoid0[1].lat+trapezoid0[2].lat+trapezoid0[3].lat)/4, lon:(trapezoid0[0].lon+trapezoid0[1].lon+trapezoid0[2].lon+trapezoid0[3].lon)/4}
let trapezoid = [];
for( let coordinate of trapezoid0){
let distanceKm = distance( coordinate.lat, coordinate.lon, gravity.lat, gravity.lon) / 1000;
let lat = coordinate.lat - (coordinate.lat - gravity.lat) * resolutionKm / 10 / distanceKm;
let lon = coordinate.lon - (coordinate.lon - gravity.lon) * resolutionKm / 10 / distanceKm;
trapezoid.push( { lat:lat, lon:lon}); // 必要以上に点が消去されることを防ぐため、予め各点を解像度(resolutionKm)の10分の1だけ重心に近づける
}
if( isInnerTriangle( point, [trapezoid[0],trapezoid[1],trapezoid[2]]) || isInnerTriangle( point, [trapezoid[0],trapezoid[2],trapezoid[3]])){
return true;
}else{
return false;
}
}
// 特定の点が三角形内に存在するかを判定
function isInnerTriangle( point, triangle){
const xr = point.lon - triangle[0].lon, yr = point.lat - triangle[0].lat;
const xp = triangle[1].lon - triangle[0].lon, yp = triangle[1].lat - triangle[0].lat;
const xq = triangle[2].lon - triangle[0].lon, yq = triangle[2].lat - triangle[0].lat;
const a = ( yr * xq - xr * yq) / ( xq * yp - xp * yq);
const b = ( xr * yp - yr * xp) / ( xq * yp - xp * yq);
if( 0<a && 0<b && a+b<1){
return true;
}else{
return false;
}
}
最後に、必須ではないかもしれませんが、ぶつ切りになった線分をなるべく結合します
// 与えられた座標の配列のうち、他の配列の先頭末尾に近いものを結合する
function jointCoordinates( coordinates, resolutionKm=5){
while( detectNearestPoint(coordinates,resolutionKm).distance < resolutionKm){
let indexes = detectNearestPoint(coordinates,resolutionKm).indexes;
let i = indexes[0], j = indexes[1], iIndex = indexes[2], jIndex = indexes[3];
if( iIndex==0){ // 1つ目の配列の先頭要素が2つ目の配列の先頭末尾に近い場合には1つ目の配列を反転(常に1つ目の配列の末尾に2つ目の配列の先頭を結合すれば良いようにする)
coordinates[i].reverse();
}
if( jIndex==-1){ // 2つ目の配列の末尾素が1つ目の配列の先頭末尾に近い場合には2つ目の配列を反転(常に1つ目の配列の末尾に2つ目の配列の先頭を結合すれば良いようにする)
coordinates[j].reverse();
}
coordinates[i] = coordinates[i].concat(coordinates[j]); // 1つ目の配列に2つ目の配列を結合
coordinates.splice( j, 1); // 結合されたほうの配列を除去
}
for( let coordinate of coordinates){
if( distance( coordinate[0].lat, coordinate[0].lon, coordinate[coordinate.length-1].lat, coordinate[coordinate.length-1].lat)/1000 < resolutionKm * 2){
coordinate.push( coordinate[0]); // 配列の先頭末尾の点が近い場合には末尾に先頭の要素を追加することで1周するようにする
}
}
return coordinates;
}
// 複数の座標の配列の先頭末尾同士で、最も近いものの距離(km)とインデックスを返す
function detectNearestPoint( coordinates0,resolutionKm=5){
let nearestDistance = Infinity, nearestIndexes = [ null, null, null, null];
// 中略。記事末尾のサンプルコードを参照
return { "distance":nearestDistance, "indexes":nearestIndexes}
}
ちゃんと計算できたか
以上の簡略化……というより大幅なサボりにより、1500行のFORTRANコードが300行のJavaScriptコードとなりました。そもそもほとんどの処理を回転楕円体近似すらせず直交直線座標系で行ったりしていますが、果たしてちゃんと書けるのでしょうか
適当な事例で描画してみると、
↓円と接線の描画まで行ったところ
↓内部の点の除去、線分の結合まで行ったところ
割ときれいに暴風警戒域を包絡する線を抽出できたように見えます
ただ、寄りで見ると一部不要な線が見えます。半径の大きな円どうしが接しているような部分では距離の変化が緩やかになってしまい、距離による内包判定が甘くなってしまう部分があるようです
というわけで、課題はありますがある程度高解像度で計算すれば簡易的な描画には使えそうです
テストページ全体のソースコード
正味の計算は calcEggplantCoordinates() 以降からです
<!DOCTYPE html>
<html lang="ja">
<head>
<meta charset="utf-8">
<title>台風(地図)</title>
<!-- https://leafletjs.com/ -->
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.9.4/dist/leaflet.css"
integrity="sha256-p4NxAoJBhIIN+hmNHrzRCf9tD/miZyoHS5obTRR9BMY="
crossorigin=""/>
<script src="https://unpkg.com/leaflet@1.9.4/dist/leaflet.js"
integrity="sha256-20nQCchB9co0qIjJZRGuk2/Z9VM+kNiyxNV1lvTlZBo="
crossorigin=""></script>
<style>
*{ font-family:sans-serif;}
html, body{ width:100%; height:100%; margin:0; overflow:hidden;}
body{ display:flex; flex-direction:column;}
div#menu{ flex:0.01;}
div#map{ flex:0.99; background:#3b4580;}
.center{ text-align:center;}
.icon{ background:#0000;}
.plain{ font-weight:bold; font-size:14px; -webkit-text-stroke:2px black; paint-order: stroke fill;}
</style>
</head>
<body>
<div id="menu">
<select id="targetTc"></select>
</div>
<div id="map"></div>
<script>
"use strict";
let layers = {'fcst':L.layerGroup()};
let map = L.map('map').setView([29, 128], 5);
initialize();
// 地図の初期化
function initialize(){
L.control.layers({
"気象庁地図":L.tileLayer('https://www.jma.go.jp/tile/jma/gray-cities/{z}/{x}/{y}.png', { minZoom:4, maxZoom:12, maxNativeZoom:10, attribution: '<a href="https://www.jma.go.jp/jma/kishou/info/coment.html">気象庁</a>'}).addTo(map),
"地理院地図":L.tileLayer('https://maps.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png', { minZoom:2, maxZoom:18, attribution: '<a href="https://maps.gsi.go.jp/development/ichiran.html">国土地理院</a>'})
}).addTo(map);
display();
}
// 表示
function display( forecast){
if( map.hasLayer(layers['fcst'])){
map.removeLayer(layers['fcst']);
}
layers['fcst'] = L.layerGroup();
let probabilityCircles = [[25.8,136.8,185],[26.0,134.3,250],[26.3,132.0,290],[26.8,128.0,320],[26.9,125.1,370],[27.5,123.1,370],[29.4,120.9,390]];
let probabilityEggplants = calcEggplantCoordinates( probabilityCircles, 5);
let probabilityEggplantLine = L.polyline(probabilityEggplants, {opacity:0, color:"#ff2800", weight:2, opacity:1, fillOpacity:0});
layers['fcst'].addLayer( probabilityEggplantLine);
map.addLayer( layers['fcst']);
}
// 複数の円の中心と半径から、それらを包摂する線の座標群を求める関数
// 引数 circles = [[lat1, lon1, r(km)], [lat2, lon2, r(km)],...]
// 返り値 coordinates = [[{lat,lon},{lat,lon},{lat,lon}...],[{lat,lon},{lat,lon},{lat,lon}...]]
function calcEggplantCoordinates( circles, resolutionKm=5){
const circleTrapezoidPoints = calcCircleTrapezoidPoints( circles, resolutionKm);
// return circleTrapezoidPoints.coordinates;
const eggplantCoordinates = removeInnerPoints( circles, circleTrapezoidPoints.trapezoids, circleTrapezoidPoints.coordinates, resolutionKm);
// return eggplantCoordinates;
return jointCoordinates( eggplantCoordinates, resolutionKm);
}
// 複数の円上の点と、となりあう円の接線上の点の座標を列挙する関数。circles = [[lat1, lon1, r(km)], [lat2, lon2, r(km)],...]
function calcCircleTrapezoidPoints( circles, resolutionKm=5){
let trapezoids = [], coordinates = [];
for( let i=0; i<circles.length; i++){
let lat1 = circles[i][0], lon1 = circles[i][1], r1 = circles[i][2];
coordinates.push( generateCircleCoordinates( lat1, lon1, r1, resolutionKm));
if( i!=circles.length-1){
let lat2 = circles[i+1][0], lon2 = circles[i+1][1], r2 = circles[i+1][2];
const tangents = getExternalTangents(lat1, lon1, r1, lat2, lon2, r2);
for( let coordinate of complementExternalTangents( tangents, resolutionKm)){
coordinates.push( coordinate);
}
trapezoids.push( [tangents[0][0], tangents[0][1], tangents[1][1], tangents[1][0]]);
}
}
return { coordinates, trapezoids};
}
// 座標群から円・台形の内部の点を除去
function removeInnerPoints( circles, trapezoids, coordinates0, resolutionKm=5){
let coordinates = [], coordinate = [];
for( let coordinate0 of coordinates0){
for( let point of coordinate0){
let lat = point.lat, lon = point.lon, isInner = false;
for( let circle of circles){
if( distance( lat, lon, circle[0], circle[1])/1000 < ( circle[2] - resolutionKm/10) ){ // 点と円の中心の距離が円の半径より小さい場合内部と判定するが、必要以上に点が消去されることを防ぐため解像度(resolutionKm)の10分の1だけ余裕をとる
isInner = true;
break;
}
}
for( let trapezoid of trapezoids){
if( isInnerTrapezoid( point, trapezoid, resolutionKm)){
isInner = true;
break;
}
}
if( isInner){
if( coordinate.length>0){
coordinates.push( coordinate);
}
coordinate = [];
}else{
coordinate.push( point);
}
}
if( coordinate.length>0){
coordinates.push( coordinate);
}
coordinate = [];
}
return coordinates;
}
// 特定の点が台形内に存在するかを判定
function isInnerTrapezoid( point, trapezoid0, resolutionKm=5){
const gravity = { lat:(trapezoid0[0].lat+trapezoid0[1].lat+trapezoid0[2].lat+trapezoid0[3].lat)/4, lon:(trapezoid0[0].lon+trapezoid0[1].lon+trapezoid0[2].lon+trapezoid0[3].lon)/4}
let trapezoid = [];
for( let coordinate of trapezoid0){
let distanceKm = distance( coordinate.lat, coordinate.lon, gravity.lat, gravity.lon) / 1000;
let lat = coordinate.lat - (coordinate.lat - gravity.lat) * resolutionKm / 10 / distanceKm;
let lon = coordinate.lon - (coordinate.lon - gravity.lon) * resolutionKm / 10 / distanceKm;
trapezoid.push( { lat:lat, lon:lon}); // 必要以上に点が消去されることを防ぐため、予め各点を解像度(resolutionKm)の10分の1だけ重心に近づける
}
if( isInnerTriangle( point, [trapezoid[0],trapezoid[1],trapezoid[2]]) || isInnerTriangle( point, [trapezoid[0],trapezoid[2],trapezoid[3]])){
return true;
}else{
return false;
}
}
// 特定の点が三角形内に存在するかを判定
function isInnerTriangle( point, triangle){
const xr = point.lon - triangle[0].lon, yr = point.lat - triangle[0].lat;
const xp = triangle[1].lon - triangle[0].lon, yp = triangle[1].lat - triangle[0].lat;
const xq = triangle[2].lon - triangle[0].lon, yq = triangle[2].lat - triangle[0].lat;
const a = ( yr * xq - xr * yq) / ( xq * yp - xp * yq);
const b = ( xr * yp - yr * xp) / ( xq * yp - xp * yq);
if( 0<a && 0<b && a+b<1){
return true;
}else{
return false;
}
}
// 与えられた座標の配列のうち、他の配列の先頭末尾に近いものを結合する
function jointCoordinates( coordinates, resolutionKm=5){
while( detectNearestPoint(coordinates,resolutionKm).distance < resolutionKm){
let indexes = detectNearestPoint(coordinates,resolutionKm).indexes;
let i = indexes[0], j = indexes[1], iIndex = indexes[2], jIndex = indexes[3];
if( iIndex==0){ // 1つ目の配列の先頭要素が2つ目の配列の先頭末尾に近い場合には1つ目の配列を反転(常に1つ目の配列の末尾に2つ目の配列の先頭を結合すれば良いようにする)
coordinates[i].reverse();
}
if( jIndex==-1){ // 2つ目の配列の末尾素が1つ目の配列の先頭末尾に近い場合には2つ目の配列を反転(常に1つ目の配列の末尾に2つ目の配列の先頭を結合すれば良いようにする)
coordinates[j].reverse();
}
coordinates[i] = coordinates[i].concat(coordinates[j]); // 1つ目の配列に2つ目の配列を結合
coordinates.splice( j, 1); // 結合されたほうの配列を除去
}
for( let coordinate of coordinates){
if( distance( coordinate[0].lat, coordinate[0].lon, coordinate[coordinate.length-1].lat, coordinate[coordinate.length-1].lat)/1000 < resolutionKm * 2){
coordinate.push( coordinate[0]); // 配列の先頭末尾の点が近い場合には末尾に先頭の要素を追加することで1周するようにする
}
}
return coordinates;
}
// 複数の座標の配列の先頭末尾同士で、最も近いものの距離(km)とインデックスを返す
function detectNearestPoint( coordinates0,resolutionKm=5){
let nearestDistance = Infinity, nearestIndexes = [ null, null, null, null];
for( let i=0; i<coordinates0.length; i++){
let iLastIndex = coordinates0[i].length-1;
for( let j=0; j<coordinates0.length; j++){
if( i==j){
continue;
}
let jLastIndex = coordinates0[j].length-1;
if( distance( coordinates0[i][0].lat, coordinates0[i][0].lon, coordinates0[j][0].lat, coordinates0[j][0].lon)/1000 < nearestDistance){
nearestDistance = distance( coordinates0[i][0].lat, coordinates0[i][0].lon, coordinates0[j][0].lat, coordinates0[j][0].lon)/1000;
nearestIndexes = [ i, j, 0, 0];
}
if( distance( coordinates0[i][0].lat, coordinates0[i][0].lon, coordinates0[j][jLastIndex].lat, coordinates0[j][jLastIndex].lon)/1000 < nearestDistance){
nearestDistance = distance( coordinates0[i][0].lat, coordinates0[i][0].lon, coordinates0[j][jLastIndex].lat, coordinates0[j][jLastIndex].lon)/1000;
nearestIndexes = [ i, j, 0, -1];
}
if( distance( coordinates0[i][iLastIndex].lat, coordinates0[i][iLastIndex].lon, coordinates0[j][0].lat, coordinates0[j][0].lon)/1000 < nearestDistance){
nearestDistance = distance( coordinates0[i][iLastIndex].lat, coordinates0[i][iLastIndex].lon, coordinates0[j][0].lat, coordinates0[j][0].lon)/1000;
nearestIndexes = [ i, j, -1, 0];
}
if( distance( coordinates0[i][iLastIndex].lat, coordinates0[i][iLastIndex].lon, coordinates0[j][jLastIndex].lat, coordinates0[j][jLastIndex].lon)/1000 < nearestDistance){
nearestDistance = distance( coordinates0[i][iLastIndex].lat, coordinates0[i][iLastIndex].lon, coordinates0[j][jLastIndex].lat, coordinates0[j][jLastIndex].lon)/1000;
nearestIndexes = [ i, j, -1, -1];
}
}
}
return { "distance":nearestDistance, "indexes":nearestIndexes}
}
// 緯度経度をローカル直交座標(m)へ変換
function latLonToXY(lat, lon, lat0, lon0) {
const R = 6378137.0;
const x = (lon - lon0) * Math.cos((lat0 * Math.PI) / 180) * Math.PI / 180 * R;
const y = (lat - lat0) * Math.PI / 180 * R;
return { x, y};
}
// ローカル直交座標(m)を緯度経度へ変換
function xyToLatLon(x, y, lat0, lon0) {
const R = 6378137.0;
const lat = lat0 + (y / R) * 180 / Math.PI;
const lon = lon0 + (x / (R * Math.cos(lat0 * Math.PI / 180))) * 180 / Math.PI;
return { lat, lon};
}
// 2円の共通外接線の接点を求める。返り値は [ [[lat,lon], [lat,lon]], ... ]
function getExternalTangents( lat1, lon1, r1, lat2, lon2, r2){
r1 = r1 * 1000;
r2 = r2 * 1000;
const originLat = lat1;
const originLon = lon1;
const c1 = { x: 0, y: 0 };
const c2 = latLonToXY( lat2, lon2, originLat, originLon);
const dx = c2.x - c1.x;
const dy = c2.y - c1.y;
const d = Math.hypot(dx, dy);
// 共通外接線が存在しない
if (d <= Math.abs(r1 - r2)) {
return [];
}
const ex = dx / d;
const ey = dy / d;
const cosTheta = (r1 - r2) / d;
if (Math.abs(cosTheta) > 1) {
return [];
}
const sinTheta = Math.sqrt(1 - cosTheta * cosTheta);
const result = [];
for (const sign of [1, -1]) {
const nx = ex * cosTheta - sign * ey * sinTheta;
const ny = ey * cosTheta + sign * ex * sinTheta;
const p1 = {
x: c1.x + r1 * nx,
y: c1.y + r1 * ny
};
const p2 = {
x: c2.x + r2 * nx,
y: c2.y + r2 * ny
};
result.push([
xyToLatLon( p1.x, p1.y, originLat, originLon),
xyToLatLon( p2.x, p2.y, originLat, originLon)
]);
}
return result;
}
// 共通外接線の接点間を補完する関数
function complementExternalTangents( tangents0, resolutionKm=5) {
let tangents = [];
for( let tangent0 of tangents0){
const distanceKm = distance( tangent0[0].lat, tangent0[0].lon, tangent0[1].lat, tangent0[1].lon) / 1000;
const points = Math.ceil( distanceKm / resolutionKm);
let tangent = [];
for (let i = 0; i < points; i++) {
tangent.push({
lat: tangent0[0].lat + (tangent0[1].lat - tangent0[0].lat) * i / points,
lon: tangent0[0].lon + (tangent0[1].lon - tangent0[0].lon) * i / points
})
}
tangents.push( tangent);
}
return tangents;
}
// 緯度経度から距離を求める関数。返り値はメートル単位。https://zenn.dev/bellbind/scraps/bb2cfb915430f5 を参考に作成
function distance( lat1, lon1, lat2, lon2){
const rad = Math.PI / 180, r = 6378137;
const p1 = lat1 * rad, p2 = lat2 * rad, q1 = lon1 * rad, q2 = lon2 * rad;
return r * Math.acos(Math.cos(p1) * Math.cos(p2) * Math.cos(q1 - q2) + Math.sin(p1) * Math.sin(p2));
}
// 特定の緯度経度を中心とした円上の各点の緯度・経度群を求める関数(中心緯度, 中心経度, 半径, 解像度)
function generateCircleCoordinates(centerLat, centerLon, radiusKm, resolutionKm = 10) {
const coords = [];
const points = Math.ceil( 2 * radiusKm * Math.PI / resolutionKm);
for (let i = 0; i < points; i++) {
const angle = (i / points) * 360;
coords.push(getCirclePoint(centerLat, centerLon, radiusKm, angle));
}
return coords;
}
// 特定の緯度経度から特定の距離・角度にある点の緯度経度を求める関数(緯度, 経度, 距離, 角度)
function getCirclePoint(centerLat, centerLon, radiusKm, angleDeg) {
const R = 6378.137; // 地球の半径 (km)
const radLat = (centerLat * Math.PI) / 180, radLon = (centerLon * Math.PI) / 180;
const radAngle = (angleDeg * Math.PI) / 180;
const dR = radiusKm / R; // 角度としての距離
// 緯度の計算
const lat2 = Math.asin(
Math.sin(radLat) * Math.cos(dR) +
Math.cos(radLat) * Math.sin(dR) * Math.cos(radAngle)
);
// 経度の計算
const lon2 =
radLon +
Math.atan2(
Math.sin(radAngle) * Math.sin(dR) * Math.cos(radLat),
Math.cos(dR) - Math.sin(radLat) * Math.sin(lat2)
);
return {
lat: (lat2 * 180) / Math.PI,
lon: (lon2 * 180) / Math.PI,
};
}
</script>
</body>
</html>