📝 Pythonで実装されたジオイド補正ツール 🐍

このツールは、国土地理院が提供するジオイド計算ツールと同じ精度と機能を持ちながら、Pythonで実装されたバージョンです。
測量やGIS（地理情報システム）を利用する専門家から、地理データの精度を向上させたい研究者まで、誰でも手軽に使えるよう設計されています。

💡 ツールの特徴
📈 ツールの使い方
📜 ソースコードの解説
実行例

💡 ツールの特徴

高精度: 国土地理院の公式ツールと同等のジオイドモデル（GSIGEO2011）を利用し、正確なジオイド補正値を提供します。
Python実装: 簡単にカスタマイズ可能で、Python環境でのスクリプト実行が可能です。
リアルタイム計算: 双一次補間法を使用し、瞬時に補正値を算出します。

📈 ツールの使い方

Python環境でスクリプトを実行します。
ジオイドデータファイル（例: gsigeo2011_ver2_2.asc）を指定します。
緯度と経度を入力して、結果としてその地点のジオイド補正値を表示します。

📜 ソースコードの解説

以下は、このツールの主要な機能を実装したPythonスクリプトです。公式ツールと同じデータを使用しながら、ユーザーが理解しやすいようにコードを整理しています。
このスクリプトは、ジオイドデータの読み込み、緯度・経度の格子点の計算、そして双一次補間法による補正値の計算という3つの主要な部分で構成されています。

1. ジオイドデータの読み込み

スクリプトは、国土地理院のデータファイルを読み込み、NumPy配列に格納します。

import numpy as np

def load_geoid_data(file_path):
    data = []
    with open(file_path, 'r') as file:
        next(file)
        buffer = []
        for line in file:
            row = [float(value) for value in line.strip().split()]
            buffer.extend(row)
            if len(buffer) == 1201:
                data.append(buffer)
                buffer = []
    if len(data) != 1801:
        raise ValueError("データの行数が不正です")
    return np.array(data)

2. 緯度と経度の格子点の計算

緯度と経度の格子点を計算し、配列として返します。

def calculate_lat_lon():
    latitudes = np.array([20.0 + (i - 1) * (1.0 / 60.0) for i in range(1, 1802)])
    longitudes = np.array([120.0 + (j - 1) * (1.5 / 60.0) for j in range(1, 1202)])
    return latitudes, longitudes

3. 補正値の計算

双一次補間法で指定した緯度と経度の補正値を計算します。

def interpolate_geoid_value(lat, lon, geoid_data, latitudes, longitudes):
    if lat < 20 or lat > 50 or lon < 120 or lon > 150:
        return 999.0000
    i = np.searchsorted(latitudes, lat) - 1
    j = np.searchsorted(longitudes, lon) - 1
    Z11 = geoid_data[i, j]
    Z12 = geoid_data[i, j + 1]
    Z21 = geoid_data[i + 1, j]
    Z22 = geoid_data[i + 1, j + 1]
    if 999.0000 in [Z11, Z12, Z21, Z22]:
        return 999.0000
    t = (lat - latitudes[i]) / (latitudes[i + 1] - latitudes[i])
    u = (lon - longitudes[j]) / (longitudes[j + 1] - longitudes[j])
    Z = (1 - t) * (1 - u) * Z11 + (1 - t) * u * Z12 + t * (1 - u) * Z21 + t * u * Z22
    return round(Z, 4)

実行例

ユーザーが緯度と経度を入力し、補正値を計算するコード例です。

input_coordinates = input("緯度と経度をカンマ区切りで入力してください（例: 35.0, 135.0）: ")
try:
    lat_str, lon_str = input_coordinates.split(',')
    lat = float(lat_str.strip())
    lon = float(lon_str.strip())
    geoid_data = load_geoid_data('gsigeo2011_ver2_2.asc')
    latitudes, longitudes = calculate_lat_lon()
    geoid_value = interpolate_geoid_value(lat, lon, geoid_data, latitudes, longitudes)
    print(f"緯度 {lat}, 経度 {lon} のジオイド補正値: {geoid_value}")
except ValueError:
    print("入力形式が正しくありません。緯度と経度はカンマ区切りで数値を入力してください。")