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

「測量」という言葉の由来 〜 測天量地 GISアプリ
「測量」という言葉の由来 〜 測天量地

📝 c#で実装されたジオイド高補正ツール

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

💡 ツールの特徴

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

📈 ツールの使い方

  1. Visual Studioや任意のC#開発環境でプロジェクトをセットアップします。
  2. ジオイドデータファイル(例: gsigeo2011_ver2_2.asc)を指定します。
  3. 緯度と経度を入力して、結果としてその地点のジオイド補正値を表示します。

📜 ソースコードの解説

1. ジオイドデータの読み込みと格子点の計算

このプログラムは、GeoidCalculatorクラスでジオイドデータの読み込みと緯度・経度の格子点の計算を行います。

以下のコードは、ジオイドデータファイルを読み込み、二次元配列に格納します。データファイルの各行を読み込んで、必要なフォーマットに変換しています。

using System;
using System.IO;
using System.Globalization;
using System.Collections.Generic;

namespace GeoidCorrection
{
    public class GeoidCalculator
    {
        private const int Nla = 1801; // 緯度方向のデータ数
        private const int Nlo = 1201; // 経度方向のデータ数
        private const double GlamN = 20.0; // 南端の緯度
        private const double GlomN = 120.0; // 西端の経度
        private const double Dgla = 1.0 / 60.0; // 緯度の間隔
        private const double Dglo = 1.5 / 60.0; // 経度の間隔
        private const double InvalidValue = 999.0000; // 無効値

        private double[,] geoData = new double[Nla, Nlo];

        // ジオイドデータをファイルから読み込む
        public bool LoadGeoidData(string filePath)
        {
            try
            {
                using (StreamReader sr = new StreamReader(filePath))
                {
                    // 1行目を読み飛ばす(ヘッダー行として扱う)
                    string headerLine = sr.ReadLine();

                    int latIndex = 0; // 緯度方向のインデックス

                    // 全データを格納するためのリスト
                    List allDataLines = new List();

                    // ファイルからすべてのデータを読み込む
                    while (!sr.EndOfStream)
                    {
                        string line = sr.ReadLine();
                        if (!string.IsNullOrWhiteSpace(line))
                        {
                            allDataLines.Add(line);
                        }
                    }

                    // 各緯度行に対して43行ずつ取り出して結合し、経度データを構成する
                    for (int i = 0; i < Nla; i++)
                    {
                        string combinedLine = string.Join(" ", allDataLines.GetRange(i * 43, 43));

                        // 行をスペースで分割
                        string[] values = combinedLine.Split(new[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
                        if (values.Length != Nlo)
                        {
                            Console.WriteLine($"データのフォーマットが不正です。行 {i + 1} に {Nlo} 個の値がありません。");
                            return false;
                        }

                        for (int lonIndex = 0; lonIndex < Nlo; lonIndex++)
                        {
                            if (!double.TryParse(values[lonIndex], NumberStyles.Float, CultureInfo.InvariantCulture, out double value))
                            {
                                Console.WriteLine($"データの変換に失敗しました: {values[lonIndex]}");
                                return false;
                            }
                            geoData[latIndex, lonIndex] = value;
                        }

                        latIndex++;
                    }
                }
                return true;
            }
            catch (Exception ex)
            {
                Console.WriteLine($"エラーが発生しました: {ex.Message}");
                return false;
            }
        }

        // ジオイド補正値を内挿計算する
        public double InterpolateGeoidValue(double latitude, double longitude)
        {
            // 緯度と経度のインデックスを計算
            int i = (int)((latitude - GlamN) / Dgla);
            int j = (int)((longitude - GlomN) / Dglo);

            if (i < 0 || i >= Nla - 1 || j < 0 || j >= Nlo - 1)
            {
                return InvalidValue; // 範囲外の場合
            }

            // 四隅の格子点の値を取得
            double z00 = geoData[i, j];
            double z01 = geoData[i, j + 1];
            double z10 = geoData[i + 1, j];
            double z11 = geoData[i + 1, j + 1];

            if (z00 == InvalidValue || z01 == InvalidValue || z10 == InvalidValue || z11 == InvalidValue)
            {
                return InvalidValue; // 無効値が含まれる場合
            }

            // tとuの計算
            double phi_i = GlamN + i * Dgla;
            double lambda_j = GlomN + j * Dglo;
            double t = (latitude - phi_i) / Dgla;
            double u = (longitude - lambda_j) / Dglo;

            // 双一次補間計算
            double z = (1 - t) * (1 - u) * z00 + (1 - t) * u * z01 + t * (1 - u) * z10 + t * u * z11;
            return Math.Round(z, 4);
        }
    }
}

このメソッドでは、データの整合性を確認し、不正なフォーマットがあればエラーメッセージを表示して処理を中断します。また、データファイルの読み込みが成功すると、内部の二次元配列に値を格納します。

さらに、格子点(緯度と経度)を計算するロジックも含まれており、格子点の間隔(緯度は1分間隔、経度は1.5分間隔)に基づいて配列を形成します。これにより、入力された緯度・経度に対して適切なデータポイントを特定し、双一次補間法を用いて補正値を算出する準備を整えます。

2. 補正値の計算

GeoidCalculatorクラスには、双一次補間法を使用して、指定した緯度と経度に基づいてジオイド補正値を計算するメソッドが含まれています。以下のコードでは、入力された座標が有効な範囲内にあるか確認し、必要に応じて補正値を返します。

using System;
using System.Globalization;
using System.IO;

namespace GeoidCorrection
{
    class Program
    {
        static void Main(string[] args)
        {
            // GeoidCalculatorのインスタンスを作成
            GeoidCalculator geoidCalculator = new GeoidCalculator();

            // データファイルのパスを指定
            string appDirectory = AppDomain.CurrentDomain.BaseDirectory;
            string filePath = Path.Combine(appDirectory, "gsigeo2011_ver2_2.asc");

            // ファイルが存在するか確認
            if (!File.Exists(filePath))
            {
                Console.WriteLine("ジオイドデータファイルが見つかりません。ファイルが正しい場所に存在するか確認してください。");
                return;
            }

            // ジオイドデータをロード
            if (!geoidCalculator.LoadGeoidData(filePath))
            {
                Console.WriteLine("ジオイドデータの読み込みに失敗しました。");
                return;
            }

            while (true)
            {
                try
                {
                    // 緯度と経度の入力を求める
                    Console.Write("緯度を入力してください: ");
                    string latText = Console.ReadLine();
                    Console.Write("経度を入力してください: ");
                    string lonText = Console.ReadLine();

                    // 緯度と経度をパース
                    double latitude = double.Parse(latText, CultureInfo.InvariantCulture);
                    double longitude = double.Parse(lonText, CultureInfo.InvariantCulture);

                    // ジオイド補正値を計算
                    double geoidValue = geoidCalculator.InterpolateGeoidValue(latitude, longitude);

                    // 結果を表示
                    if (geoidValue == 999.0000)
                    {
                        Console.WriteLine("指定された座標のジオイド補正値は無効です。");
                    }
                    else
                    {
                        Console.WriteLine($"緯度 {latitude}, 経度 {longitude} のジオイド補正値: {geoidValue:F4} m");
                    }

                    // 再度入力を求めるかどうか
                    Console.WriteLine("続けて入力しますか?(y/n): ");
                    string continueInput = Console.ReadLine();
                    if (continueInput.ToLower() != "y")
                    {
                        break;
                    }
                }
                catch (FormatException)
                {
                    Console.WriteLine("入力形式が正しくありません。数値を入力してください。");
                }
                catch (Exception ex)
                {
                    Console.WriteLine($"エラーが発生しました: {ex.Message}");
                }
            }
        }
    }
}

このメソッドでは、以下のような手順で計算を行います:

  • 入力された緯度・経度が有効範囲に収まっているか確認します。
  • 格子点に基づき、周辺4点のデータを取得します。
  • 双一次補間法により、指定された座標での補正値を算出します。

実行例

ユーザーが緯度と経度を入力し、補正値を計算するコード例です。実行時にコンソールから緯度と経度を入力し、リアルタイムでジオイド補正値を表示します。

    // コンソールアプリケーション形式で、ユーザーに対して緯度・経度の入力を促します。
    // その後、入力された値をもとにジオイド補正値を計算し、結果を表示します。
    


コメント

タイトルとURLをコピーしました