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

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

ツールの特徴

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

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

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

ツールの使い方

Visual Studioや任意のC++開発環境でプロジェクトをセットアップします。

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

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

ソースコードの解説

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

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

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

#include <cstdio>

#include <cmath>

#include <cstring>



const int MAX_ROW = 1801;  // 正しい行数

const int MAX_COL = 1201;  // 正しい列数

const double ERROR_VALUE = 999.0;

const double EL2 = 0.00001;

const double DX = 1.5 / 60.0;

const double DY = 1.0 / 60.0;

const double X_MIN = 120.0;

const double Y_MIN = 20.0;



// グローバル変数

static double gHeight[MAX_ROW][MAX_COL] = { { ERROR_VALUE } };



// ファイルからジオイド高を読み込む関数

extern "C" __declspec(dllexport) int __stdcall Geoid_Read(LPCTSTR fin) {

    FILE* fp = fopen(fin, "r");

    if (fp == nullptr) {

        return VB_FALSE;

    }



    char str[254];

    double wHeight[29] = { 0 };

    int row = 0;



    // ファイルを行ごとに読み込む

    while (fgets(str, sizeof(str), fp) != nullptr) {

        // 1行目はスキップ

        if (row < 1) {

            row++;

            continue;

        }



        // ジオイド高をパース

        for (int L = 1; L <= 28; L++) {

            int start = (L - 1) * 9;

            char buf[10] = { 0 };

            strncpy(buf, str + start, 9);

            wHeight[L] = atof(buf);

        }



        int I = (row - 1) / 43 + 1;

        int L = (row - 1) % 43;

        if (L == 0) {

            L = 43;

            I--;

        }



        int J = (L - 1) * 28;

        for (int L = 1; L <= 28; L++) {

            gHeight[I][L + J] = wHeight[L];

        }



        row++;

    }



    fclose(fp);

    return 0;

}



// 線形補間を用いてジオイド高を求める関数

extern "C" __declspec(dllexport) double __stdcall Geoid_Bilinear(double lat, double lon) {

    double xPt = lon;

    double yPt = lat;



    long ix = static_cast<long>((xPt - X_MIN) / DX) + 1;

    long iy = static_cast<long>((yPt - Y_MIN) / DY) + 1;



    double x = (xPt - X_MIN) / DX - (ix - 1);

    double y = (yPt - Y_MIN) / DY - (iy - 1);



    long jx = ix + 1;

    long jy = iy + 1;



    // データ範囲外のチェック

    if (ix < 0 || ix >= MAX_COL || iy < 0 || iy >= MAX_ROW) {

        return ERROR_VALUE;

    }



    // グリッド上かをチェック

    int iadx = (std::fabs(x) < EL2) ? 0 : ((1.0 - std::fabs(x)) < EL2) ? 1 : 99;

    int iady = (std::fabs(y) < EL2) ? 0 : ((1.0 - std::fabs(y)) < EL2) ? 1 : 99;



    if (iady < 10) {

        if (iadx < 10) {

            return gHeight[iy + iady][ix + iadx];

        } else {

            if (gHeight[iy + iady][ix] == ERROR_VALUE || gHeight[iy + iady][jx] == ERROR_VALUE) {

                return ERROR_VALUE;

            }

            return (1.0 - x) * gHeight[iy + iady][ix] + x * gHeight[iy + iady][jx];

        }

    } else if (iadx < 10) {

        if (gHeight[iy][ix + iadx] == ERROR_VALUE || gHeight[jy][ix + iadx] == ERROR_VALUE) {

            return ERROR_VALUE;

        }

        return (1.0 - y) * gHeight[iy][ix + iadx] + y * gHeight[jy][ix + iadx];

    }



    // グリッド線上にない点の場合の補間

    if (gHeight[jy][ix] == ERROR_VALUE || gHeight[jy][jx] == ERROR_VALUE || gHeight[iy][ix] == ERROR_VALUE || gHeight[iy][jx] == ERROR_VALUE) {

        return ERROR_VALUE;

    }



    return (1.0 - x) * (1.0 - y) * gHeight[iy][ix] +

           y * (1.0 - x) * gHeight[jy][ix] +

           x * (1.0 - y) * gHeight[iy][jx] +

           x * y * gHeight[jy][jx];

}