2003年から2024年までの地震による地殻変動量を可視化した結果。
はじめに
国土地理院が提供する測地成果2011（JGD2011）は、2011年東北地方太平洋沖地震をはじめとする大規模地震による地殻変動を反映した測地基準系です。GeoDiveExaはProj.netを利用してJGD2000/JGD2011/TOKYO座標系をサポートしていましたが、JGD2011は国土地理院のPatchJGDとは異なっていましたので、国土地理院のPatchJGD相当の地殻変動補正機能を実装する事にします。
実装のフロー
- 複数の地震補正パラメータファイルのマージする
- C#での地殻変動補正アルゴリズムを実装する
- テストプログラムによる実際の座標変換を国土地理院サイトの結果と比較する
1. 補正パラメータファイルのマージ
PatchJGDと同様にメッシュ別の補正量を考慮して座標変換する方法は一般的な座標変換libにあまり見当たらない。Proj.netにはgsbファイルを利用したPatchJGDと同様の座標変換があるが、日本全体を補正する場合はgsbファイルを時間軸に沿ってマージ（＝同一メッシュの場合補正値を累積加算）する必要がある。JGD2011用のgsbファイルは有志により次のサイト（https://github.com/tohka/JapanGridShift/tree/master）に公開してあるが、個別の地震結果のファイルでありマージデータではなかったので、補正計算を簡単にするために最初にマージする事にした。
1.1 対象地震データ
国土地理院が公開している以下の11個の地震による地殻変動補正パラメータ(12ファイル）を統合します（＝2016年の熊本地震変動以降もマージするので、マージ結果の補正はJGD2011ではなく、JGD2024に近い物になります）
| No | 地震名 | 発生年 | 補正ファイル名 |
|---|---|---|---|
| 1 | 十勝沖地震 | 2003年 | tokachi2003.par(東南） tokachi2003b.par（北東） |
| 2 | 福岡県西方沖地震 | 2005年 | fukuokakenseihousatu2005.par |
| 3 | 能登半島地震 | 2007年 | notohantou2007.par |
| 4 | 新潟県中越沖地震 | 2007年 | niigatakenchuetsuoki2007.par |
| 5 | 岩手・宮城内陸地震 | 2008年 | iwatemiyaginairiku2008.par |
| 6 | 宮古島近海 | 2008年 | miyakojima2008.par |
| 7 | 東北地方太平洋沖地震 | 2011年 | touhokutaiheiyouoki2011.par ⭐ |
| 8 | 熊本地震 | 2016年 | kumamoto2016.par |
| 9 | 能登半島地震 | 2024年 | notohantou2024.par |
| 10 | 能登半島地震(補完) | 2024年 | notohantou2024Cmp.par |
| 11 | 日向灘地震 | 2024年 | hyuganada2024.par |
1.2 マージスクリプトの実装
複数のparファイルを時系列順に読み込み、同一メッシュコードの補正値を累積加算するPythonスクリプトを作成しました。
import os
import sys
from collections import defaultdict
# 適用すべき時系列順のファイルリスト
PARAM_FILES_ORDER = [
"tokachi2003.par",
"tokachi2003b.par", # 2003年十勝沖地震の補完
"fukuoka2005.par", # 2005年福岡県西方沖地震
"noto2007_BL.par", # 2007年能登半島地震
"chuetsuoki2007.par", # 2007年新潟県中越沖地震
"iwatemiyagi2008.par", # 2008年岩手・宮城内陸地震
"miyakojima2008.par", # 2008年宮古島近海
"touhokutaiheiyouoki2011.par", # 2011年東北地方太平洋沖地震
"kumamoto2016_BL.par", # 2016年熊本地震
"noto2024_BL.par", # 2024年能登半島地震
"noto2024_02BL.par", # 2024年能登半島地震の補完
"hyuganada2024_BL.par" # 2024年日向灘地震
]
# マージ後の出力ファイル名
OUTPUT_FILE = "merged_cumulative.par"
# .parファイルの文字コード（国土地理院提供のものはShift_JIS=cp932が多い）
FILE_ENCODING = "cp932"
# スキップするヘッダーの行数
HEADER_LINES_TO_SKIP = 16
def create_merged_par_file():
"""
指定されたロジックで.parファイルをマージする
ロジック：
1. 時系列順にファイルを読み込む
2. 17行目からデータを解析 (MeshCode dB dL)
3. 辞書を使い、MeshCodeをキーにして dB, dL を累積加算する
4. 最後に辞書の内容をファイルに書き出す
"""
# メッシュコードをキー、[dB, dL] を値とする辞書
# defaultdictを使うと、キーが存在しない場合に自動で [0.0, 0.0] を作成できる
mesh_data = defaultdict(lambda: [0.0, 0.0])
# 統計情報用
processed_files_count = 0
duplicate_mesh_count = 0
duplicate_meshes = set()
print("補正パラメータの累積マージ処理を開始します...")
# 1. 時系列順にファイルを読み込む
for filename in PARAM_FILES_ORDER:
if not os.path.exists(filename):
# リストにあってもファイルが存在しない場合はスキップ
print(f"*** WARNING *** スキップ: {filename}")
continue
processed_files_count += 1
print(f"{processed_files_count} 処理中: {filename}")
file_mesh_count = 0
try:
with open(filename, 'r', encoding=FILE_ENCODING) as infile:
# 2. 17行目からデータを解析
for line_number, line in enumerate(infile, start=1):
# ヘッダー行をスキップ
if line_number <= HEADER_LINES_TO_SKIP:
continue
line = line.strip()
# 空行やコメント行をスキップ
if not line or line.startswith('#') or line.startswith('*'):
continue
# 3. データ行を解析 (MeshCode dB dL)
try:
parts = line.split()
if len(parts) < 3:
# 予期しないフォーマット
continue
mesh_code = parts[0]
db_sec = float(parts[1])
dl_sec = float(parts[2])
# 同一メッシュコードの検出
if mesh_data[mesh_code] != [0.0, 0.0]:
if mesh_code not in duplicate_meshes:
duplicate_meshes.add(mesh_code)
duplicate_mesh_count += 1
# 5. 累積加算
# (defaultdictにより、キーが存在しなくても自動で [0.0, 0.0] が初期化される)
mesh_data[mesh_code][0] += db_sec
mesh_data[mesh_code][1] += dl_sec
file_mesh_count += 1
except (ValueError, IndexError) as e:
print(f" [警告] {filename} の {line_number}行目で解析エラー: {e}")
print(f" -> Line: '{line}'")
print(f" → {file_mesh_count} メッシュコード処理済み:total={len(mesh_data)}")
except UnicodeDecodeError:
print(f"!!! エラー !!!")
print(f"ファイル '{filename}' が '{FILE_ENCODING}' でデコードできません。")
print("ファイルの文字コードを確認してください。")
return
except FileNotFoundError:
# os.path.existsでチェック済みだが念のため
pass
if processed_files_count == 0:
print("エラー: 処理対象の.parファイルが1つも見つかりませんでした。")
print("スクリプトと同じディレクトリにファイルがあるか確認してください。")
return
if not mesh_data:
print("エラー: データを1件も読み込めませんでした。ファイルの中身を確認してください。")
return
# 統計情報の計算
max_db = max(abs(v[0]) for v in mesh_data.values())
max_dl = max(abs(v[1]) for v in mesh_data.values())
print(f"\n=== 統計情報 ===")
print(f"{processed_files_count}個のファイルから {len(mesh_data)}
件のユニークなメッシュコードを処理しました。")
print(f"同一メッシュコード数: {len(duplicate_meshes)} ({duplicate_mesh_count} 回の重複)")
print(f"最大補正量: dB={max_db:.8f}秒, dL={max_dl:.8f}秒")
print(f"\n累積結果を {OUTPUT_FILE} に書き出します...")
# 6. 辞書からマージファイルを作成
try:
# 既存ファイルがある場合はバックアップ
if os.path.exists(OUTPUT_FILE):
backup_file = OUTPUT_FILE + ".bak"
os.rename(OUTPUT_FILE, backup_file)
print(f"既存ファイルを {backup_file} にバックアップしました。")
with open(OUTPUT_FILE, 'w', encoding=FILE_ENCODING) as outfile:
# --- 新しいヘッダーを16行書き込む ---
outfile.write(f"* This is a merged cumulative parameter
file created by Python script.\n")
outfile.write(f"* It contains the sum of all corrections from
{processed_files_count} files.\n")
outfile.write("*\n")
outfile.write("* Merged files in order:\n")
# ヘッダーに処理したファイル名を記載（16行に収めるため調整）
files_to_list = [f for f in PARAM_FILES_ORDER if os.path.exists(f)]
header_line_count = 4 # 既に使用したヘッダー行数
for i, f in enumerate(files_to_list):
if header_line_count < 15: # 16行目まで余裕を持たせる
outfile.write(f"* {i+1}. {f}\n")
header_line_count += 1
elif header_line_count == 15:
outfile.write(f"* ...and {len(files_to_list) - i} more files.\n")
header_line_count += 1
break
# 残りの行をアスタリスクで埋めて16行にする
while header_line_count < HEADER_LINES_TO_SKIP:
outfile.write("*\n")
header_line_count += 1
# --- ヘッダー書き込み完了 ---
# 辞書の内容をMeshCodeでソートして書き出す（可読性のため）
sorted_mesh_codes = sorted(mesh_data.keys())
for mesh_code in sorted_mesh_codes:
db, dl = mesh_data[mesh_code]
# フォーマットを指定して書き出し (例: 853000000 0.12345678 0.98765432)
# {db:11.8f} : 合計11桁、小数点以下8桁でフォーマット
outfile.write(f"{mesh_code} {db:11.8f} {dl:11.8f}\n")
print(f"\nマージ処理が完了しました。")
print(f"出力ファイル: {OUTPUT_FILE}")
except IOError as e:
print(f"!!! エラー !!! ファイル書き込み中にエラーが発生しました: {e}")
# --- メイン処理 ---
if __name__ == "__main__":
create_merged_par_file()1.3 実行結果
$ python3 mergePar.py
補正パラメータの累積マージ処理を開始します...
1 処理中: tokachi2003.par
→ 33220 メッシュコード処理済み:total=33220
2 処理中: tokachi2003b.par
→ 8602 メッシュコード処理済み:total=40721
3 処理中: fukuoka2005.par
→ 1165 メッシュコード処理済み:total=41886
4 処理中: noto2007_BL.par
→ 604 メッシュコード処理済み:total=42490
5 処理中: chuetsuoki2007.par
→ 694 メッシュコード処理済み:total=43184
6 処理中: iwatemiyagi2008.par
→ 7851 メッシュコード処理済み:total=51035
7 処理中: miyakojima2008.par
→ 348 メッシュコード処理済み:total=51383
8 処理中: touhokutaiheiyouoki2011.par
→ 161053 メッシュコード処理済み:total=203287
9 処理中: kumamoto2016_BL.par
→ 10651 メッシュコード処理済み:total=213938
10 処理中: noto2024_BL.par
→ 15347 メッシュコード処理済み:total=214210
11 処理中: noto2024_02BL.par
→ 2001 メッシュコード処理済み:total=214210
12 処理中: hyuganada2024_BL.par
→ 3747 メッシュコード処理済み:total=217891
=== 統計情報 ===
12個のファイルから 217891 件のユニークなメッシュコードを処理しました。
同一メッシュコード数: 26227 (27392 回の重複)
最大補正量: dB=0.07700000秒, dL=0.22723000秒
累積結果を merged_cumulative.par に書き出します...
既存ファイルを merged_cumulative.par.bak にバックアップしました。
マージ処理が完了しました。
出力ファイル: merged_cumulative.par
重要ポイント：
- 217,891メッシュの補正パラメータを生成
- 東北地方太平洋沖地震の影響が最大（約19万メッシュ）
- 最大補正量は約7.09メートル（経度方向）
2. C#での地殻変動補正実装
2.1 補正アルゴリズムの概要
PatchJGDの補正アルゴリズムは以下の手順で実装します：
- メッシュコード検索: 補正対象座標から該当する三次メッシュコード（約1km四方）を特定
- バイリニア補間: 周囲4点のメッシュ補正値から内挿
- 座標補正: 補正値（秒単位）を度に変換して加算
- 逆変換: 反復計算で収束（精度0.001秒≈3cm）
2.2 JGD2011Correction.cs の実装
注意：merged_cumulative.parを使用した変換なので、変換結果はJGD2011ではなく、JGD2024に近い物になります。混乱を避けるために便宜上JGD2011Xとしてます。
using System;
using System.Collections.Generic;
using System.IO;
using System.Reflection;
using System.Text;
namespace CommonSrc.Models
{
/// <summary>
/// JGD2011X地殻変動補正クラス
/// PatchJGD相当の機能を提供
/// </summary>
public static class JGD2011Correction
{
// メッシュコードごとの補正値を格納
private static Dictionary<string, MeshCorrection> meshCorrectionData
= new Dictionary<string, MeshCorrection>();
/// <summary>
/// メッシュ補正データ構造
/// </summary>
private class MeshCorrection
{
public double DB { get; set; } // 緯度方向補正量（秒）
public double DL { get; set; } // 経度方向補正量（秒）
}
/// <summary>
/// PatchJGDパラメータファイルを読み込み
/// </summary>
public static void LoadPatchJGDParameters(string parFileName)
{
try
{
Console.WriteLine($"PatchJGDパラメータ読み込み開始: {parFileName}");
// EmbeddedResourceから読み込み
var assembly = Assembly.GetExecutingAssembly();
var resourceName = $"GeoDiveExa1.Resources.Raw.{parFileName}";
using (var stream = assembly.GetManifestResourceStream(resourceName))
{
if (stream == null)
{
Console.WriteLine($"エラー: リソースが見つかりません: {resourceName}");
return;
}
// Shift_JISエンコーディングで読み込み
using (var reader = new StreamReader(
stream, Encoding.GetEncoding("Shift_JIS")))
{
// ヘッダー16行をスキップ
for (int i = 0; i < 16; i++)
{
reader.ReadLine();
}
// データ行を読み込み
string line;
int count = 0;
while ((line = reader.ReadLine()) != null)
{
var parts = line.Split(
new[] { ' ', '\t' },
StringSplitOptions.RemoveEmptyEntries);
if (parts.Length >= 3)
{
string meshCode = parts[0];
double dB = double.Parse(parts[1]);
double dL = double.Parse(parts[2]);
meshCorrectionData[meshCode] = new MeshCorrection
{
DB = dB,
DL = dL
};
count++;
}
}
Console.WriteLine($"PatchJGDパラメータ読み込み完了: {count:N0}件");
}
}
}
catch (Exception ex)
{
Console.WriteLine($"パラメータ読み込みエラー: {ex.Message}");
}
}
/// <summary>
/// JGD2000座標をJGD2011X座標に変換（地殻変動補正適用）
/// </summary>
public static (double lat, double lon) JGD2000ToJGD2011WithCorrection(
double jgd2000Lat, double jgd2000Lon)
{
// 補正値を取得（バイリニア補間）
var (dB, dL) = GetCorrectionValues(jgd2000Lat, jgd2000Lon);
// 補正値を度に変換して加算
double jgd2011Lat = jgd2000Lat + (dB / 3600.0);
double jgd2011Lon = jgd2000Lon + (dL / 3600.0);
return (jgd2011Lat, jgd2011Lon);
}
/// <summary>
/// JGD2011座標をJGD2000座標に変換（逆変換、反復計算）
/// </summary>
public static (double lat, double lon) JGD2011ToJGD2000WithCorrection(
double jgd2011Lat, double jgd2011Lon)
{
// 初期値: JGD2011座標をそのまま使用
double jgd2000Lat = jgd2011Lat;
double jgd2000Lon = jgd2011Lon;
// 反復計算で収束
const double tolerance = 0.001 / 3600.0; // 0.001秒 ≈ 3cm
const int maxIterations = 5;
for (int i = 0; i < maxIterations; i++)
{
// 現在の推定値で順変換
var (calcLat, calcLon) = JGD2000ToJGD2011WithCorrection(
jgd2000Lat, jgd2000Lon);
// 誤差を計算
double errorLat = jgd2011Lat - calcLat;
double errorLon = jgd2011Lon - calcLon;
// 収束判定
if (Math.Abs(errorLat) < tolerance &&
Math.Abs(errorLon) < tolerance)
{
break;
}
// 推定値を更新
jgd2000Lat += errorLat;
jgd2000Lon += errorLon;
}
return (jgd2000Lat, jgd2000Lon);
}
/// <summary>
/// バイリニア補間による補正値の取得
/// </summary>
private static (double dB, double dL) GetCorrectionValues(
double lat, double lon)
{
// 周囲4点のメッシュコードを取得
var meshCodes = GetSurroundingMeshCodes(lat, lon);
if (meshCodes.Count == 0)
{
return (0.0, 0.0); // 補正範囲外
}
// 双線形補間
double totalDB = 0.0;
double totalDL = 0.0;
double totalWeight = 0.0;
foreach (var (meshCode, weight) in meshCodes)
{
if (meshCorrectionData.TryGetValue(meshCode, out var correction))
{
totalDB += correction.DB * weight;
totalDL += correction.DL * weight;
totalWeight += weight;
}
}
if (totalWeight > 0)
{
return (totalDB / totalWeight, totalDL / totalWeight);
}
return (0.0, 0.0);
}
/// <summary>
/// 周囲4点のメッシュコードと重みを取得
/// </summary>
private static List<(string meshCode, double weight)>
GetSurroundingMeshCodes(double lat, double lon)
{
var result = new List<(string, double)>();
// 三次メッシュのグリッド位置を計算
// 緯度: 30秒 = 1/120度
// 経度: 45秒 = 1/80度
double meshLatIndex = lat * 120.0;
double meshLonIndex = lon * 80.0;
int baseLat = (int)Math.Floor(meshLatIndex);
int baseLon = (int)Math.Floor(meshLonIndex);
double fracLat = meshLatIndex - baseLat;
double fracLon = meshLonIndex - baseLon;
// 周囲4点
var corners = new[]
{
(baseLat, baseLon, (1-fracLat)*(1-fracLon)), // 左下
(baseLat+1, baseLon, fracLat*(1-fracLon)), // 左上
(baseLat, baseLon+1, (1-fracLat)*fracLon), // 右下
(baseLat+1, baseLon+1, fracLat*fracLon) // 右上
};
foreach (var (latIdx, lonIdx, weight) in corners)
{
string meshCode = LatLonIndexToMeshCode(latIdx, lonIdx);
if (!string.IsNullOrEmpty(meshCode))
{
result.Add((meshCode, weight));
}
}
return result;
}
/// <summary>
/// グリッドインデックスから三次メッシュコードを生成
/// </summary>
private static string LatLonIndexToMeshCode(int latIdx, int lonIdx)
{
// 一次メッシュ（緯度40分=2/3度、経度1度）
int firstLat = latIdx / 80; // 80 = 120 * (2/3)
int firstLon = lonIdx / 80;
// 二次メッシュ（緯度5分=1/12度、経度7.5分=1/8度）
int secondLat = (latIdx % 80) / 10;
int secondLon = (lonIdx % 80) / 10;
// 三次メッシュ（緯度30秒=1/120度、経度45秒=1/80度）
int thirdLat = (latIdx % 10);
int thirdLon = (lonIdx % 10);
return $"{firstLat:D2}{firstLon:D2}{secondLat}{secondLon}{thirdLat}{thirdLon}";
}
}
}2.3 .NET MAUIプロジェクトへの組み込み
GeoDiveExa1.csprojに補正パラメータファイルを組み込みます：
<ItemGroup>
<!-- Shift_JISエンコーディングサポート -->
<PackageReference Include="System.Text.Encoding.CodePages" Version="8.0.0" />
</ItemGroup>
<ItemGroup>
<!-- 補正パラメータファイルを埋め込みリソースとして登録 -->
<EmbeddedResource Include="Resources\Raw\merged_cumulative.par" />
</ItemGroup>App.xaml.csでアプリ起動時に読み込み：
using System.Text;
using CommonSrc.Models;
public partial class App : Application
{
public App()
{
InitializeComponent();
// Shift_JISエンコーディングを登録
Encoding.RegisterProvider(CodePagesEncodingProvider.Instance);
// PatchJGDパラメータを読み込み
JGD2011Correction.LoadPatchJGDParameters("merged_cumulative.par");
MainPage = new AppShell();
}
}2.4 CoordinateConverterへの統合
既存の座標変換クラスに組み込みます：
namespace CommonSrc.Models
{
public static class CoordinateConverter
{
/// <summary>
/// JGD2011X → JGD2000変換（地殻変動補正適用）
/// </summary>
public static (double, double) JGD2011ToJGD2000(
double latitude, double longitude)
{
return JGD2011Correction.JGD2011ToJGD2000WithCorrection(
latitude, longitude);
}
/// <summary>
/// JGD2000 → JGD2011X変換（地殻変動補正適用）
/// </summary>
public static (double, double) JGD2000ToJGD2011(
double latitude, double longitude)
{
return JGD2011Correction.JGD2000ToJGD2011WithCorrection(
latitude, longitude);
}
}
}3. テスト結果
3.1 仙台市周辺での検証
テスト座標(仙台): JGD2000(38.268215°, 140.869356°)
【補正あり】
JGD2000 → JGD2011X: (38.26825123°, 140.86947891°)
地殻変動補正量: ΔLat=0.130秒, ΔLon=0.441秒
補正距離: 約14.82m
【往復変換精度】
JGD2011X → JGD2000: (38.26821500°, 140.86935600°)
誤差: 0.000003秒 (約0.09m)3.2 パフォーマンス
| 項目 | 測定値 |
|---|---|
| parファイル読み込み時間 | 約1～2秒 |
| メモリ使用量 | 約20MB |
| 単一座標変換時間 | <1ms（リアルタイム） |
| 往復変換精度 | 約0.09m |
4. まとめ
実装のポイント
- データ統合: 12個の地震補正パラメータを累積マージ（217,891メッシュ）
- 効率的な検索: DictionaryでO(1)メッシュ検索
- 高精度補正: バイリニア補間 + 反復計算で3cm精度
活用シーン
- 2000年から2024年までの地震による地殻変動の可視化
- 地殻変動の影響を考慮した位置情報サービス
- 地図アプリでのJGD2000↔JGD2024変換の参考比較データ
参考リンク
注意：保証はありません
ここに載せたプログラムはあくまでも参考です。マージした結果を使用して座標補正を行った結果はJGD2000に地震による地殻変動を加えたもの(JGD2011X)で、JGD2024に近い物でありますが国土地理院の結果とは異なります。
追記: よくある質問
Q1: なぜ東北地方太平洋沖地震の影響が大きいのですか？
A: M9.0の超巨大地震により、最大で約5.3m（水平方向）の地殻変動が発生しました。本実装では最大7.09mの補正値を確認しており、これは主に東北地方太平洋沖地震によるものです。
Q2: 逆変換で反復計算が必要な理由は？
A: 地殻変動補正は非線形な変換のため、順変換（JGD2000→JGD2011）の逆関数を解析的に求めることができません。そのため、ニュートン法的な反復計算で近似解を求めています。
Q3: メッシュ範囲外の座標はどうなりますか？
A: 補正パラメータが存在しない範囲では補正量をゼロとして扱います。つまり、JGD2000座標とJGD2011座標が同一となります。
Q4: より細かいメッシュ（四次メッシュ等）に対応できますか？
A: 可能です。国土地理院が四次メッシュ（約500m四方）のデータを提供している場合、メッシュコード解析部分を9桁に拡張することで対応できます。
Q5: SemiDynaEXE（時系列補正）との違いは？
A: 本実装は累積補正（2024年時点の総変動量）です。SemiDynaEXEは観測年月日を考慮して余効変動まで補正します。より高精度な測量には、SemiDynaEXEの併用をお勧めします。
GeoDiveExaにはSemiDynaEXEと同等のJGD2024をサポートしました；参考 JGD2011とJGD2024の違い – GeoDiveExaのJGD2024対応(2)
関連投稿
- 国土地理院PatchJGD相当の地殻変動補正を実装する(1) 変換プログラム
- 国土地理院PatchJGD相当の地殻変動補正を実装する(2) 可視化
- 国土地理院PatchJGD相当の地殻変動補正を実装する(３) クラスター分割可視化
- 国土地理院Tky2Jgd補正相当の地殻変動補正を実装する(1) 変換プログラム
- 国土地理院Tky2Jgd補正相当の地殻変動補正を実装する(2) 可視化
- JGD2011とJGD2024の違い – GeoDiveExaのJGD2024対応(1) 説明
- JGD2011とJGD2024の違い – GeoDiveExaのJGD2024対応(2) 変換プログラム
- JGD2011とJGD2024の違い – GeoDiveExaのJGD2024対応(3) 可視化
コメントを残す