Zenn
📑

走時表をもとに時間と震源の深さからP波とS波の大きさを求めてみよう(地図に円を描くところまで)

2023/04/19に公開
C#
Leaflet
#
地震
tech

この記事は 地震界隈 Advent Calendar 2020 10日目の記事です。
"ただただ参加したかった" というだけなのでとんだクソ記事ですが許してください。
技術力がないので大した記事が書けませんが、ふと思い浮かんだものを書きたいと思います。

なにをするの？

緊急地震速報で発表される震源の深さと発生時刻と現在時刻の差を用いてP波/S波の予報円の大きさを求めます。（タイトルそのまま）
なお、線形補間を使用しているため誤差等がありますがそこは目を瞑ってください。あくまでも趣味の範疇です。（逃げ）

後半ではサンプルとして Web のマップシステムとしてメジャーな Leaflet を使用して実際に地図に描いてみようと思います。

おことわり

なお、筆者は数学やその他の知識が乏しく誤りなどがあるかもしれません。あった場合はコメント等で "やさしく" 教えてください。

走時表って？

これ です。
走時と距離の関係を表にしたものですね。はい。

本題

本セクション中のコードは C# です。

走時表準備

ただのテキストファイルなので簡単に取り込めます。改行区切りで配列にしておきましょう。

Program.cs
private class TravelTimeTable
{
    public double P { get; set; }
    public double S { get; set; }
    public int Depth { get; set; }
    public int Distance { get; set; }
}

このようなクラスを用意しておきます。

Program.cs
private static TravelTimeTable[] _table;

public static async void ImportData() {
    var file = await File.ReadAllTextAsync("tjma2001"); // ファイルを読み込む
    _table = Regex.Replace(file, "\x20+", " ") // 連続した空白を一つの空白に置換
      .Trim() // 前後の空白を除去
      .Replace("\r", "") // \r\n を \n へ
      .Split('\n') // 改行ごとに分割
      .Select(x => x.Split(' ')) // それぞれを空白ごとに分割
      .Select(x =>
          new TravelTimeTable // 型を変換しそれぞれを格納
          {
              P = double.Parse(x[1]),
              S = double.Parse(x[3]),
              Depth = int.Parse(x[4]),
              Distance = int.Parse(x[5])
          })
      .ToArray(); // TravelTimeTable[] へ変換
}

これで _table に走時表の各行のデータが格納されました。

円の大きさを求める

では実際に求めてみましょう。
流れとしては

  1. 渡された深さのデータに絞る(この時点でレコードが残らなかったら終了)
  2. 渡された時間をもとに それ以上の最初のデータ と それ以下の最後のデータ を抽出する
  3. 線形補間で大きさを求める

となっています。これをP波とS波の両方で行います。

  • int depth は深さを km で
  • int time は発生時刻と求めたい時刻の差を秒で
Program.cs
public static (double, double) GetValue(int depth, int time)
{
  // 走時表では "深さ701km以深" "2001秒以上" は計算できないため
  if (depth > 700 || time > 2000) return (double.NaN, double.NaN);

  var values = _table.Where(x => x.Depth == depth).ToArray(); // 1
  if (values.Length == 0) return (double.NaN, double.NaN);

  var p1 = values.LastOrDefault(x => x.P <= time); // 2
  var p2 = values.FirstOrDefault(x => x.P >= time); // 2
  if (p1 == null || p2 == null) return (double.NaN, double.NaN);
  var p = (time - p1.P) / (p2.P - p1.P) * (p2.Distance - p1.Distance) + p1.Distance; // 3

  var s1 = values.LastOrDefault(x => x.S <= time); // 2
  var s2 = values.FirstOrDefault(x => x.S >= time); // 2
  if (s1 == null || s2 == null) return (p, double.NaN);
  var s = (time - s1.S) / (s2.S - s1.S) * (s2.Distance - s1.Distance) + s1.Distance; // 3

  return (p, s);
}

Tuple で返していますが、 item1 がP波の半径(km)で item2 がS波の半径(km)です。

動作確認

動作確認環境では TravelTimeTableConverter というクラスにまとめました。

まず最初に走時表を読み込んでから計算させます。

Program.cs
private static void Main()
{
  TravelTimeTableConverter.ImportData();
  Console.WriteLine("distance time : (p(km) , s(km))");
  Console.WriteLine($"20km 20sec : {TravelTimeTableConverter.GetValue(20, 20)}");
  Console.WriteLine($"100km 200sec : {TravelTimeTableConverter.GetValue(100, 200)}");
  Console.WriteLine($"200km 200sec : {TravelTimeTableConverter.GetValue(200, 200)}");
  Console.WriteLine($"300km 200sec : {TravelTimeTableConverter.GetValue(300, 200)}");
}
Result
$ dotnet run
distance time : (p(km) , s(km))
20km 20sec : (122.35900962861072, 67.68853695324285)
100km 200sec : (1603.2552083333333, 868.2417083144026)
200km 200sec : (1639.8745519713261, 874.7576045627376)
300km 200sec : (1672.7323943661972, 869.2659627953747)

しっかりと結果が出力されれば OK です。

これを作るきっかけをくれた 某つよいプログラマー様 が作られた 変換プログラム と比較してみても

Result
$ dotnet run
(122.35900962861074, 67.68853695324283)
(1603.2552083333333, 868.2417083144026)
(1639.8745519713261, 874.7576045627377)
(1672.7323943661972, 869.2659627953745)

おそらく正しく計算できていると思われます。（自信なし）

短いし寂しいのでサンプルをつくる

実際に使う機会があるのは地図に予報円を描画する際だと思います。それ以外は思いつきませんが。
拙作のツールの中に 緊急地震速報一覧でデータを選択したら地図に発表当時の予報円を描く というものがありますのでそれを使いたいと思います。

※画像には緊急地震速報一覧がありますが本サンプルは地図に円を描くだけです。

使用するライブラリは jQuery axios Leaflet です。

コードはこちらに置いてあります。
iedred7584/eew-map-sample-with-travel-time-table-converter

index.html
<!DOCTYPE html>
<html lang="ja">

<head>
  <meta charset="UTF-8">
  <title>eew-map-sample-with-travel-time-table-converter</title>
  <link rel="stylesheet" href="https://unpkg.com/leaflet@1.7.1/dist/leaflet.css" />
  <link rel="stylesheet" href="css/style.min.css">
</head>

<body>
  <div id="map"></div>
  <script src="https://code.jquery.com/jquery-3.5.1.min.js"></script>
  <script src="https://cdn.jsdelivr.net/npm/axios/dist/axios.min.js"></script>
  <script src="https://unpkg.com/leaflet@1.7.1/dist/leaflet.js"></script>
  <script async src="scripts/main.js"></script>
</body>

</html>
style.scss
body {
  margin: 0;
  width: 100%;
  height: 100vh;
  overflow: hidden;
}

#map {
  width: 100%;
  height: 100vh;
}

JavaScript へメソッド名を変えずに移植して計算させます。

main.js
(async () => {
  // 地図を用意
  const map = new L.Map("map", {
    center: new L.LatLng(38.558595, 137.0550225),
    zoom: 6,
    preferCanvas: true
  });

  // 適当なタイルマップを設定
  L.tileLayer("https://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png").addTo(map);

  // 走時表を読み込む
  const table = await ImportTable();

  // 計算
  const value = GetValue(table, 60, 57);

  //L.circle の第二引数は Radius (半径) を m で指定するため * 1000 をする
  // P波
  L.circle([35.6896342, 139.6921007], value[0] * 1000, {
    color: "#0000ff",
    weight: 3,
    opacity: 1,
    fillColor: '#000000',
    fillOpacity: 0
  }).addTo(map);

  // S波
  L.circle([35.6896342, 139.6921007], value[1] * 1000, {
    color: "#ff0000",
    weight: 6,
    opacity: 1,
    fillColor: '#ff0000',
    fillOpacity: 0.2
  }).addTo(map);

  // 震央もどき
  L.marker([35.6896342, 139.6921007], {
    icon: L.icon({
      iconUrl: "assets/epicenter.png",
      iconSize: [60, 60],
      iconAnchor: [30, 30]
    })
  }).addTo(map)
})()

async function ImportTable() {
  return (await axios.get("../assets/tjma2001.txt")).data
    .trim().replace("\r", "").replace(/\x20+/g, " ").split("\n")
    .map(x => {
      const s = x.split(" ");
      return {
        p: parseFloat(s[1]),
        s: parseFloat(s[3]),
        depth: parseInt(s[4]),
        distance: parseInt(s[5]),
      };
    });
}

function GetValue(table, depth, time) {
  if (depth > 700 || time > 2000) return [NaN, NaN];

  const values = table.filter(x => x.depth == depth);
  if (values.Length == 0) return [NaN, NaN];

  let p1 = values.filter(x => x.p <= time);
  p1 = p1[p1.length - 1];
  let p2 = values.filter(x => x.p >= time);
  p2 = p2[0];
  if (p1 == null || p2 == null) return [NaN, NaN];
  const p = (time - p1.p) / (p2.p - p1.p) * (p2.distance - p1.distance) + p1.distance;

  let s1 = values.filter(x => x.s <= time);
  s1 = s1[s1.length - 1];
  let s2 = values.filter(x => x.s >= time);
  s2 = s2[0];
  if (s1 == null || s2 == null) return [p, NaN];
  const s = (time - s1.s) / (s2.s - s1.s) * (s2.distance - s1.distance) + s1.distance;

  return [p, s];
}

完成です。

Discussion