国土地理院PatchJGD相当の地殻変動補正を実装する(2)

前回マージして作成した日本全体に対する地殻変動補正パラメータファイルが正しいのか？確認するために変動量を画像表示してみます。

2003年から2024年までの地震による地殻変動量を可視化した結果。

１補正データの可視化

1.1 可視化スクリプトの作成

217,891メッシュの補正データを可視化するPythonスクリプトを作成しました。

主な機能：

ヒートマップ表示（緯度・経度・総移動量）
移動ベクトル図
クラスタリングと外接線表示
ゼロ値の特別表示（黒色）

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
JGD2011補正パラメータの可視化スクリプト
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
from sklearn.cluster import DBSCAN

# 設定
INPUT_FILE = "merged_cumulative.par"
OUTPUT_FOLDER = "image"
CLUSTER_EPS = 0.3  # 約33km
CLUSTER_MIN_SAMPLES = 50

def create_clusters_and_hulls(data):
    """DBSCANクラスタリングと凸包計算"""
    coords = data[:, [1, 0]]  # (lon, lat)

    clustering = DBSCAN(
        eps=CLUSTER_EPS, 
        min_samples=CLUSTER_MIN_SAMPLES
    ).fit(coords)

    labels = clustering.labels_
    hulls = []

    for label in set(labels):
        if label == -1:  # ノイズをスキップ
            continue

        cluster_coords = coords[labels == label]
        if len(cluster_coords) >= 4:
            try:
                hull = ConvexHull(cluster_coords)
                hulls.append(cluster_coords[hull.vertices])
            except:
                pass

    return labels, hulls

def plot_heatmap_with_clusters(data, values, title, filename):
    """クラスター境界線付きヒートマップ"""
    fig, ax = plt.subplots(figsize=(12, 10))

    # 散布図でメッシュデータを表示
    scatter = ax.scatter(
        data[:, 1], data[:, 0],  # lon, lat
        c=values,
        s=1,
        marker='s',
        cmap='RdBu_r'
    )

    plt.colorbar(scatter, ax=ax, label='移動量 (秒)')

    # クラスター境界線を描画
    labels, hulls = create_clusters_and_hulls(data)
    for hull_points in hulls:
        hull_closed = np.vstack([hull_points, hull_points[0]])
        ax.plot(
            hull_closed[:, 0], hull_closed[:, 1],
            color='darkblue',
            linewidth=2.0,
            alpha=0.7
        )

    ax.set_xlabel('経度 (度)')
    ax.set_ylabel('緯度 (度)')
    ax.set_title(f'{title}\nクラスター数: {len(hulls)}')
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')

    plt.savefig(f'{OUTPUT_FOLDER}/{filename}', dpi=150, bbox_inches='tight')
    plt.close()

    print(f'保存完了: {filename} (クラスター数: {len(hulls)})')