• 
• トップ ›    日記/一般｜  浜松市｜　
• 新規登録｜ログイン
• はまぞう検索⇒

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

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

ツールの特徴

• 高精度: 国土地理院の公式ツールと同等のジオイドモデル（GSIGEO2011）を利用し、正確なジオイド補正値を提供します。


• Python実装: 簡単にカスタマイズ可能で、Python環境でのスクリプト実行が可能です。


• リアルタイム計算: 双一次補間法を使用し、瞬時に補正値を算出します。

ツールの使い方

1. Python環境でスクリプトを実行します。


2. ジオイドデータファイル（例: gsigeo2011_ver2_2.asc）を指定します。


3. 緯度と経度を入力して、結果としてその地点のジオイド補正値を表示します。

ソースコードの解説

以下は、このツールの主要な機能を実装した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("入力形式が正しくありません。緯度と経度はカンマ区切りで数値を入力してください。")