TKY2JGDとは。
キーワードをたどってこのページを読まれている方はご存知の方がほとんどだとは思いますが、国土地理院が配布する日本測地系という座標系の緯度・経度を、現在業界標準で使われている世界測地系という緯度・経度の座標系に変換する機能を提供するソフトウェアです。
https://www.gsi.go.jp/sokuchikijun/tky2jgd.html
元々フリーウェアとして上記サイト上で配布されていたのですが、元のソースがVBで1990年代のソフトウェアということですでにかなり古くなっており、脆弱性が見つかったことでバイナリの配布はすでに終了し、現在はソースコードの配布と、Web版の利用のみ続けられている状況です。
https://vldb.gsi.go.jp/sokuchi/surveycalc/tky2jgd/main.html
なぜC#に？
私がかかわっているプロジェクトで上記ソフトウェアを利用しているのですが、Web版でも一括変換機能が提供されているものの1度に2万件までという制限があり、配布されていた実行ファイルを利用して処理を行っていました。
ところが、このローカル版の配布が終了。また配布されているソースのプロジェクトがVB5.0で開発されていたものということで実行環境も古くなっており、Windows10では動作しないという状況になりました。
ご存知の通り、先日Window7のサポートが終了したことで、社内マシンをWindows10環境に入れ替える必要が出てきたことから、Windows10対応のツールを作成する必要が出てきたため、配布されているソースをもとにC#に置き換える対応を実施したものです。
同種のことでお困りの方がいるかなとGitHubでの公開も考えたのですが
- おそらく需要は日本国内だけ
- 社内向けツールのため（下の前提条件にも書いていますが）公開対象は動作するプロジェクト全体ではなく、TKY2JGDから置き換えた1ファイルのみ
と考えたことから、こちらに張っておけば、必要な方には届くだろうと判断し、こちらにソースコードを貼ることで公開としました。
ご利用の際には適宜、ご自身のプロジェクト内に取り込む形で修正・改修を行ってご利用ください。
なお、自分が利用する範囲内で計算結果が従来のソフトウェアと変わらないことの確認は行っていますが、計算結果の正当性や、本コードを利用して発生する不利益に対するいかなる補償もいたしかねますので、ご利用の際は自己責任においてご利用いただきますよう、ご承知ください。
前提条件
今回、私の業務に沿った対応のみ行っており、以下のような形で置き換えを行っています。
- 私は元々、座標空間系の処理を専門とする者ではないため、座標計算の正当性については判断がついていません。このため、座標処理の部分についてはVBの計算式をC#に置き換えるということだけに注力しました。（結果が、旧来のソフトウェアの変換結果と同じ値になることをもって正と判断しています）
- 修正対象としたのは、元のプロジェクトの中で座標変換のコア計算を行っている modTky2jgd.bas のみです。
- 各メソッド冒頭には、C#の標準（と思われる）形式に書き直した形で元々記載されていたコメントを記述しています。また[]内に、VBで利用されていた元のメソッド名を記載していますので、処理内容を確認する際にご利用ください。
- 処理ロジック部分については、例えば
文字列を走査して文字数をカウントしながら分割処理を行っているような箇所は Splitメソッドに置き換えるなど、できるだけわかりやすい形への置き換えを行ったつもりです。 - 各部分のコメントについても処理内容の正当性が判断つかないため、（一部手直しはしているもののの）基本的には元ソースにあったものをそのまま残しています。
- 変数名、メソッド名のコーディングルールはある程度VBからC#向けのものに置き換えました。ただ、そもそもC#もあまり得意ではないため、あまりC#らしからぬコーディングになっている箇所が多々あるかと思いますがお目こぼしいただければと思います。
- 元のソースでは、ソフトウェアのフォームオブジェクトにテキストボックスが配置され、各種数値が入力されていたことから、その入力値が適正かどうか判断するロジックが入っていましたが、これらの処理は不要（あらかじめチェック済みの前提）としています。
- 元のソースでは、画面上でいくつか設定を行い、その設定値に合わせて動的に入力・出力を判断していましたが、今回、私のプロジェクトで必要だったのが
緯度経度ファイル一括出力ファイルに入力値を出力するの設定だったため、 これらのロジック以外の部分は検証は行っておりません 。（パラメータを定数として設定しているため、Visual Studioでソースを読み込むと到達できないコードがみつかるのもそのためです）。 - 変換に必要な地域パラメータファイルは、国土地理院の配布サイトから取得してください。このファイルの現状では、全国版の
TKY2JGD.parファイルを、同一プログラム内のリソースとして保存してあるものを読み込むようにしていますが、特定地域のみにしたり、またプロジェクト内に保存せずに動的リソースとして読み込むような場合は、コンストラクタ部分とreadParameters()メソッド内の処理を適宜書き換えてお使いください。 - ソース内の
MainForm mainFormCommonUtilsConstString当たりは、私が作成しているプロジェクト部分になります。画面表示やログ出力、メッセージ用の文字列定数等々となりますので、適宜書き換えてご利用ください。
ということで、置き換え後のファイルがこちらになります。
Tky2jgd.cs
/***
* Original Program TKY2JGD https://www.gsi.go.jp/sokuchikijun/tky2jgd_download.html
* witten by (C) Mikio TOBITA 飛田幹男，国土地理院
* with Visual Basic 5.0
*
* convert to C# authering by Tezuka Yasuto 手塚泰斗
* 2020-03-24
*
***/
using MeiboConverter;
using System;
using System.IO;
using System.Reflection;
using System.Text;
using System.Windows.Forms;
namespace Tky2jgd
{
class Tky2jgd
{
private MainForm mainForm;
private readonly string inFileName; // 入力ファイル名，出力ファイル名
private readonly string outFileName; // 入力ファイル名，出力ファイル名
private readonly string parFileName; // 変換パラメータファイルのフルパス名
private readonly string parVersion; // 返還パラメータファイルのVersion番号
// パラメータファイル読み込み
private long[] paramMeshCode; //1行1ｸﾞﾘｯﾄﾞ型ﾃﾞｰﾀ変換ﾊﾟﾗﾒｰﾀを格納
private float[] paramDiffBeta; //1行1ｸﾞﾘｯﾄﾞ型ﾃﾞｰﾀ変換ﾊﾟﾗﾒｰﾀを格納
private float[] paramDiffLambda; //1行1ｸﾞﾘｯﾄﾞ型ﾃﾞｰﾀ変換ﾊﾟﾗﾒｰﾀを格納
private struct MeshCode
{
public int FirstMeshCode { get; set; } //１次メッシュコード 4桁 ex.5438
public int SecondMeshCode { get; set; } //２次メッシュコード 2桁 ex.23
public int ThridMeshCode { get; set; } //３次メッシュコード 2桁 ex.45
public long FullMeshCode { get; set; } //すべて ex.54382345
public string FormedMeshCode { get; set; } //正式な形式の文字列型 ex."5438-23-45"
public double ModLat { get; set; } //余り。3次メッシュの左下（南西角）点からどれくらいずれているか。0(西端)以上1(東端)未満
public double ModLon { get; set; } //余り。3次メッシュの左下（南西角）点からどれくらいずれているか。0(南端)以上1(北端)未満
}
private struct TranPara //座標変換パラメータ緯度差，経度差の構造体を定義
{
public double DiffBeta { get; set; }
public double DiffLambda { get; set; }
}
private static readonly double PI = 3.14159265358979; // 円周率
private static readonly double DoublePI = 6.28318530717959; // 円周率*2
private static readonly double HalfPI = 1.5707963267949; // 円周率/2
private static readonly double RadToDeg = 57.2957795130823; //ラジアンを度にするため掛ける数
private static readonly double DegToRad = 1.74532925199433E-02; //度をラジアンにするため掛ける数
private static readonly double ROBYO = 206264.806247; //3600*180/PI
private struct EllipParam //楕円体パラメータの構造体を定義
{
public double Axix { get; set; } //semimajor axix in meter
public double Flattening1 { get; set; } //reciprocal flattening
public double Flattening { get; set; } //flattening
public double Eccentricity { get; set; } //squared eccentricity
public string Namec { get; set; } //ellipsoid name
}
private EllipParam[] ellipObj; // EP[0]はBessel楕円体=日本測地系，EP[1]はGRS-80楕円体=世界測地系
private static readonly bool allThreeParamFlag = false; //すべて3parameterによる計算を行うかどうか。今回はfalse固定
private static readonly bool inputOutFlag = true; //出力ファイルに入力情報を記録するかのフラグ。元々は画面で切り替えできるが今回は true 固定
private static readonly bool inputOptR3Flag = true; //設定画面の設定 「地域ごとのパラメータ」が無かったら、「3パラメータ」で変換する。今回はfalse固定
/// <summary>
/// コンストラクタ
/// </summary>
/// <param name="targetDir">入力ファイル、出力ファイルのフォルダパス</param>
/// <param name="thisForm">表示している画面フォーム。各種情報を表示する際の操作に使用</param>
/// <param name="inFile">入力ファイル。デフォルトは既存プログラムに従いzahyou.in</param>
/// <param name="outFile">出力ファイル。デフォルトは既存プログラムに従いzahyou.out</param>
/// <param name="parmFile">パラメータファイル。デフォルトは既存プログラムの全国版でTKY2JGD.par。リソースから読み込んで処理する場合は指定は不要</param>
public Tky2jgd( string targetDir, MainForm thisForm, string inFile = "\\zahyou.in", string outFile= "\\zahyou.out", string parmFile= "\\TKY2JGD.par")
{
// ファイル名関連の固定値設定
inFileName = targetDir+ inFile;
outFileName = targetDir+ outFile;
parFileName = Application.StartupPath + parmFile;
//楕円体ﾊﾟﾗﾒｰﾀの設定
//楕円体パラメータ。0：ベッセル楕円体，1：GRS-80楕円体
ellipObj = new EllipParam[2];
ellipObj[0].Axix = 6377397.155;
ellipObj[0].Flattening1 = 299.152813;
ellipObj[0].Namec = "Bessel";
ellipObj[0].Flattening = 1.0 / ellipObj[0].Flattening1;
ellipObj[0].Eccentricity = (2.0 * ellipObj[0].Flattening1 - 1.0) / ellipObj[0].Flattening1 / ellipObj[0].Flattening1; //=e*e (squared e);
ellipObj[1].Axix = 6378137.0;
ellipObj[1].Flattening1 = 298.257222101;
ellipObj[1].Namec = "GRS-80";
ellipObj[1].Flattening = 1.0 / ellipObj[1].Flattening1;
ellipObj[1].Eccentricity = (2.0 * ellipObj[1].Flattening1 - 1.0) / ellipObj[1].Flattening1 / ellipObj[1].Flattening1; //=e*e (squared e)
mainForm = thisForm;
// パラメータファイルの読み込み
parVersion = ReadParamFile();
}
/// <summary>
/// 地域ごとパラメータファイル(TKY2JGD.par)の読み込み
/// 読み込みに成功したら""でなく，Versionを返す
/// </summary>
/// <returns></returns>
public string ReadParamFile()
{
// 1行1ｸﾞﾘｯﾄﾞ型ﾊﾟﾗﾒｰﾀﾌｧｲﾙを読み込みます。
// ﾌｧｲﾙ形式例２（２行抜粋）：
// 3次ﾒｯｼｭｺｰﾄﾞ dB(") dL(")
// 46305526 12.74179 -8.18083
// 46305624 12.75169 -8.20710
long lineCount = 0;
string firstLine = "";
string secondLine = "";
// リソースの中に保存した TKY2JGD.par ファイルを使う場合のロジック
var paramFile = MeiboConverter.Properties.Resources.TKY2JGD;
string paramLines = Encoding.GetEncoding("shift_jis").GetString(paramFile);
string[] lineArray = paramLines.Split('\n');
paramMeshCode = new long[lineArray.Length];
paramDiffBeta = new float[lineArray.Length];
paramDiffLambda = new float[lineArray.Length];
foreach (string line in lineArray)
{
switch (lineCount)
{
case 0:
firstLine = line;
break;
case 1:
secondLine = line;
break;
default:
if (line.Length == 0)
{
continue;
}
// parファイルの読み込み行。Debug用処理だが有効にしていると処理にかなり時間がかかるのでコメント
// line.DebugLog(nameof(line));
string[] param_array = line.Replace(" ", " ").Replace(" ", " ").Split(' ');
paramMeshCode[lineCount - 2] = long.Parse(param_array[0]);
paramDiffBeta[lineCount - 2] = float.Parse(param_array[1]);
paramDiffLambda[lineCount - 2] = float.Parse(param_array[2]);
if (lineCount % 10000 == 0)
{
mainForm.AddProgressValue(1);
mainForm.SetProcessMsg(string.Format(ConstString.Process.MsgReadParamFile, lineCount));
}
break;
}
lineCount += 1;
}
/** TKY2JGD.par ファイルを、どこかのフォルダから動的に読み込む場合はこちらのロジックを使う
StreamReader paramf = new StreamReader(parfilename);
try
{
FileInfo fi = new FileInfo(parfilename);
//ファイルのサイズを取得
long filesize = fi.Length;
// プログレスバーのMax値は，ﾌｧｲﾙの大きさ(in Bytes)とする。
//frmParaFile1.ProgressBar1.Max = filesize / (double)128; // for ProgressBar
// frmParaFile1.lblLine = string.Format(filesize.ToString(), "###,###,###") + " Bytes";
first_line = paramf.ReadLine();
second_line = paramf.ReadLine();
while (paramf.Peek() != -1)
{
string line = paramf.ReadLine();
line_count += 1;
string[] line_array = line.Split(' ');
MC[line_count] = long.Parse(line_array[0]);
dB1[line_count] = float.Parse(line_array[1]);
dL1[line_count] = float.Parse(line_array[2]);
if ((line_count % 2000 == 0))
{
mainForm.setProcessMsg(ConstString.Process.msgReadParamFile);
//frmParaFile1.ProgressBar1.Value = filesize;
if ((BreakFlag))
{
BreakFlag = false;
return "";
}
}
}
}
catch (Exception e)
{
CommonUtils.logWrite(string.Format(ConstString.Error.fileNotExist, parfilename));
MessageBox.Show(
ConstString.Error.fatalMsg,
string.Format(ConstString.Error.fatalTitle),
MessageBoxButtons.OK,
MessageBoxIcon.Error
);
}
finally
{
paramf.Close();
}
**/
// parファイルの1行目からVersion番号を取り出す。20～28カラムです。
return firstLine.Substring(19, 9);
}
/// <summary>
/// 順方向の座標変換をﾌｧｲﾙ一括処理します[BLFile()]
/// 座標変換後の緯度経度をdmsフォーマット（度分秒形式）で出力します。
/// </summary>
public void FileConvertBatch()
{
Assembly assembly = Assembly.GetExecutingAssembly();
AssemblyName asmName = assembly.GetName();
Version app_version = asmName.Version;
// 詳細表示をしない設定とします
//frmAdv.chkDetail.Value = 0;
// 処理中止ボタンを表示
//frmFile.cmdBreak.Visible = true;
StreamReader reader = new StreamReader(inFileName);
try
{
StreamWriter writer = new StreamWriter(outFileName, false, Encoding.GetEncoding("shift_jis"));
writer.WriteLine($"このファイル{Path.GetFileNameWithoutExtension(outFileName)}は，プログラムTKY2JGD Ver.{app_version}が{Path.GetFileNameWithoutExtension(inFileName)}を読み込んで計算処理したものです。");
writer.WriteLine("使用した変換パラメータファイルは，" + Path.GetFileNameWithoutExtension(parFileName) + " " + parVersion + "です。");
if (inputOutFlag)
{
writer.WriteLine("次に示すように，各行の最初の２つの数字が入力した日本測地系の緯度，経度，次の２つが変換されたJGD2000系の緯度，経度を表しています。");
writer.WriteLine(" 日本測地系（入力値） JGD2000系（計算値）");
writer.WriteLine(" 緯度 経度 緯度 経度");
writer.WriteLine("ddmmss.sssss dddmmss.sssss ddmmss.sssss dddmmss.sssss");
}
else
{
writer.WriteLine("次に示すように，各行の最初の２つの数字が，変換されたJGD2000系の緯度，経度を表しています。");
writer.WriteLine(" JGD2000系（計算値）");
writer.WriteLine(" 緯度 経度");
writer.WriteLine(" ddmmss.sssss dddmmss.sssss");
}
writer.WriteLine("行末に「3parameters」があるものは，地域毎のパラメータがなかったか3パラメータで変換するよう設定されていたため，「東京大正」三角点における3パラメータで変換したものです。");
writer.WriteLine("また，「-9999.」がある行は，変換されなかった行です。");
writer.WriteLine("以上のどちらでもない行は，「地域毎の変換パラメータ」で変換された行です。");
writer.WriteLine("ただし，コメント行や数値の形式が不正な行は，変換されずにそのまま出力されます。");
writer.Close();
long lineCount = 0;
while (reader.Peek() != -1) // 行ごとの処理
{
string lines = reader.ReadLine(); //行を変数に読み込みます。
#if DEBUG
CommonUtils.LogWrite("読み込みLINE ==> ["+ lines + "]");
#endif
//改行(CR/LF)だけの行または#から始まる行（コメント）行の場合、そのままの内容を出力
if (lines.Length == 0 || lines.StartsWith("#"))
{
writer.WriteLine(lines);
lineCount += 1;
continue;
}
string[] lineArray = lines.Split(' ');
// 半角スペースで分割し、要素数が２または3以外（3はコメント）の行があったらフォーマット不正で終了
if( lineArray.Length != 2 && lineArray.Length != 3)
{
CommonUtils.LogWrite( string.Format(ConstString.Error.FormatFailedLatLon, lineCount, lines) );
throw new Exception(string.Format(ConstString.Error.FormatFailedLatLon, lineCount, lines));
}
// コメント要素が無かった場合は空白を設定
if( lineArray.Length == 2)
{
Array.Resize(ref lineArray, lineArray.Length + 1);
lineArray[2] = "";
}
ConvertTkyToJgd(lineArray);
lineCount += 1;
if( lineCount%10 == 0)
{
mainForm.SetProcessMsg(string.Format(ConstString.Process.MsgFamilySplitting, lineCount));
mainForm.AddProgressValue(10);
}
}
}
catch(Exception e)
{
string Msg = e.Message;
CommonUtils.LogWrite(Msg);
MessageBox.Show(
Msg,
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Error);
throw new Exception(ConstString.Error.ErrorEndConvert);
}
finally
{
reader.Close();
}
}
/// <summary>
/// 位置データ(緯度，経度)を読み込みます。[LatLonOnly]
/// </summary>
/// <param name="inputLatLon">入力された緯度・経度の文字列</param>
private void ConvertTkyToJgd( string[] inputLatLon )
{
inputLatLon.DebugLog(nameof(inputLatLon));
string outLat = "";
string outLon = "";
double x1 = 0.0;
double y1 = 0.0;
double x2 = 0.0;
double y2 = 0.0;
int sign = 0;
string signString = "";
int degree = 0;
int minute = 0;
double second = 0.0;
// 3ﾊﾟﾗﾒｰﾀで変換したことの印
string threeParaMark = "";
double lat = Convert.ToDouble(inputLatLon[0]);
double lon = Convert.ToDouble(inputLatLon[1]);
// 計算開始
// 順方向（日本測地系からJGD2000）変換
// （１）緯度，経度を度分秒形式を度にし，
// （２）パラメータを参照し，bilinear補間法により補間し，
// （３）補間した，dL,dBを実用成果のL,Bに加える。
// （４）結果をテキストボックスに書き込む。
// （１）緯度，経度を度分秒形式を度にし，
DmsToDeg(lat, ref y1, ref sign, ref degree, ref minute, ref second);
DmsToDeg(lon, ref x1, ref sign, ref degree, ref minute, ref second);
// (A) もし，「すべて3ﾊﾟﾗﾒｰﾀで変換」なら
if (allThreeParamFlag)
{
Tokyo97toITRF94(y1, x1, ref y2, ref x2);
threeParaMark = " 3parameters";
}
else
{
// （２）パラメータを参照し，bilinear補間法により補間し，
TranPara tranPara = Bilinear1(y1, x1);
if (tranPara.DiffLambda == -9999.0 || tranPara.DiffBeta == -9999.0)
{
if (inputOptR3Flag)
{
Tokyo97toITRF94(y1, x1, ref y2, ref x2);
threeParaMark = " 3parameters";
}
else
{
threeParaMark = "";
WriteOutFile( lat, lon, outLat, outLon, inputLatLon[2], threeParaMark );
return;
}
}
else
{
// （３）補間した，dL,dBを実用成果のL,Bに加える
x2 = x1 + tranPara.DiffLambda / 3600.0;
y2 = y1 + tranPara.DiffBeta / 3600.0;
threeParaMark = "";
}
}
// （４）結果をファイルに書き込む。
DegToDms(y2, ref degree, ref minute, ref second, ref signString);
outLat = $"{signString}{degree,2}{minute:D2}{second:00.00000}";
DegToDms(x2, ref degree, ref minute, ref second, ref signString);
outLon = $"{signString}{degree,2}{minute:D2}{second:00.00000}";
#if DEBUG
CommonUtils.LogWrite($"変換後緯度・経度 ==> [{outLat},{outLon}]");
#endif
WriteOutFile( lat, lon, outLat, outLon, inputLatLon[2], threeParaMark );
}
/// <summary>
/// 結果を出力ファイルに書き出します。[Kekka]
/// </summary>
/// <param name="lat">緯度</param>
/// <param name="lon">経度</param>
/// <param name="latout">出力用にフォーマットを整えた緯度文字列</param>
/// <param name="lonout">出力用にフォーマットを整えた経度文字列</param>
/// <param name="comments">コメント</param>
/// <param name="markThreePara">3パラメータ変換したか否かのマーク文字列</param>
private void WriteOutFile(
double lat,
double lon,
string latout,
string lonout,
string comments,
string markThreePara
)
{
string out_line = "";
if (inputOutFlag)
{
// dms
string latin = $"{lat:000000.00000}";
string lonin = $"{lon:000000.00000}";
out_line = $"{latin} {lonin} {latout} {lonout} {comments}{markThreePara}\n";
}
else
{
out_line = $"{latout} {lonout} {comments}{markThreePara}\n";
}
File.AppendAllText(outFileName, out_line, Encoding.GetEncoding("shift_jis"));
}
/// <summary>
/// Bilinear interpolationをする準備とBilinear interpolationのcall
/// 1行1グリッド型変換パラメータファイル対応｡
/// 1行にdB,dLの2つの変換ﾊﾟﾗﾒｰﾀがあるので，緯度，経度，同時に補間計算を行う。
/// </summary>
/// <param name="lat">緯度 or y 単位は度</param>
/// <param name="lon">経度 or x 単位は度</param>
/// <returns></returns>
private TranPara Bilinear1(double lat, double lon)
{
MeshCode eastMC = new MeshCode();
MeshCode northMC = new MeshCode();
MeshCode northEastMC = new MeshCode();
TranPara rtnTmpPara = new TranPara();
// from Ver.1.3.78 2001/11/22
// マイナスの緯度（南緯），経度（西経）等日本の領土外は標準地域メッシュコードが定義されていない場所では，
// 地域毎の変換パラメータがないので，直ちにOutsideにとぶ。
// この判定がなかったVer.1.3.77では，メッシュコードを検索に行ってしまい。見つかってしまうバグがあった。
// なお，このバグは日本領土内では結果にまったく影響しない。
if (lat < 20 || lat > 46)
{
rtnTmpPara.DiffBeta = -9999.0;
rtnTmpPara.DiffLambda = -9999.0;
return rtnTmpPara;
}
if (lon < 120 || lon > 154)
{
rtnTmpPara.DiffBeta = -9999.0;
rtnTmpPara.DiffLambda = -9999.0;
return rtnTmpPara;
}
// lat,lonからMesh Codeを作成。
// lat,lonを含むメッシュのﾒｯｼｭｺｰﾄﾞを計算
MeshCode originMC = LatLon2MeshCode(lat, lon);
// MC0の東，北，東北隣のﾒｯｼｭｺｰﾄﾞを計算
AdjacentMeshCode( originMC, ref eastMC, ref northMC, ref northEastMC);
// parameter search。1行1ｸﾞﾘｯﾄﾞ型。データファイルの座標系つまり旧東京測地系(=旧日本測地系)で。
// 4つの内最初の1つは，ﾊﾟﾗﾒｰﾀﾌｧｲﾙの最初からサーチするが，他の3つはそれ以降順番に探すのみ。
// ﾊﾟﾗﾒｰﾀﾌｧｲﾙは低緯度から高緯度に，同じ緯度なら，低経度から高経度の順番に並ぶべき。
// メッシュコードの小さい順ではない。
// 原点検索
int search_result = 0;
search_result = Array.IndexOf(paramMeshCode, originMC.FullMeshCode);
if( search_result == -1)
{
rtnTmpPara.DiffBeta = -9999.0;
rtnTmpPara.DiffLambda = -9999.0;
return rtnTmpPara;
}
long paramSuff00 = search_result;
// 東隣
search_result = 0;
search_result = Array.IndexOf(paramMeshCode, eastMC.FullMeshCode);
if (search_result == -1)
{
rtnTmpPara.DiffBeta = -9999.0;
rtnTmpPara.DiffLambda = -9999.0;
return rtnTmpPara;
}
long paramSuff10 = search_result;
// 北隣
search_result = 0;
search_result = Array.IndexOf(paramMeshCode, northMC.FullMeshCode);
if (search_result == -1)
{
rtnTmpPara.DiffBeta = -9999.0;
rtnTmpPara.DiffLambda = -9999.0;
return rtnTmpPara;
}
long paramSuff01 = search_result;
// 北東隣
search_result = 0;
search_result = Array.IndexOf(paramMeshCode, northEastMC.FullMeshCode);
if (search_result == -1)
{
rtnTmpPara.DiffBeta = -9999.0;
rtnTmpPara.DiffLambda = -9999.0;
return rtnTmpPara;
}
long paramSuff11 = search_result;
// 変換ﾊﾟﾗﾒｰﾀの代入
// 00=原点、10=東隣、01=北隣、11=北東隣
double dB00 = paramDiffBeta[paramSuff00];
double dB10 = paramDiffBeta[paramSuff10];
double dB01 = paramDiffBeta[paramSuff01];
double dB11 = paramDiffBeta[paramSuff11];
double dL00 = paramDiffLambda[paramSuff00];
double dL10 = paramDiffLambda[paramSuff10];
double dL01 = paramDiffLambda[paramSuff01];
double dL11 = paramDiffLambda[paramSuff11];
// bilinear補間。データファイルの座標系つまり日本測地系で。
// 結果は，緯度，経度のJGD2000－日本測地系,または，Tokyo97-日本測地系
// ModLat，ModLonは0以上1未満
rtnTmpPara.DiffBeta = Interpol(dB00, dB10, dB01, dB11, originMC.ModLon, originMC.ModLat);
rtnTmpPara.DiffLambda = Interpol(dL00, dL10, dL01, dL11, originMC.ModLon, originMC.ModLat);
// grid点でのデータの表示。データファイルの座標系つまり日本測地系で。
return rtnTmpPara;
}
//
/// <summary>
/// Bilinear interpolation [interpol1]
/// X0to1, Y0to1は0以上1未満
///
/// Y↑
/// u3 u4
///
/// u1 u2 → X
/// </summary>
/// <param name="u1">原点</param>
/// <param name="u2">東隣</param>
/// <param name="u3">北隣</param>
/// <param name="u4">北東隣</param>
/// <param name="x0to1">x座標の移動距離</param>
/// <param name="y0to1">y座標の移動距離</param>
/// <returns></returns>
public double Interpol(double u1, double u2, double u3, double u4, double x0to1, double y0to1)
{
double origin;
double valueX;
double valueY;
double valueXY;
origin = u1;
valueX = u2 - u1;
valueY = u3 - u1;
valueXY = u4 - u2 - u3 + u1;
double rtn = origin + valueX * x0to1 + valueY * y0to1 + valueXY * x0to1 * y0to1;
return rtn;
}
/// <summary>
/// XYからXYへの順方向（旧日本測地系から日本測地系2000へ）の変換を行います[doXyInterpol1]
/// </summary>
/// <param name="convflag">地域ごとのパラメータで変換するか否かのフラグ</param>
/// <param name="inLat">入力緯度</param>
/// <param name="inLon">入力経度</param>
/// <param name="inX">入力x</param>
/// <param name="inY">入力y</param>
public void ForwardXYInterpol(bool convflag, double inLat, double inLon, double inX, double inY)
{
// 平面直角座標値XY（旧日本測地系）をテキストボックスから読み込み緯度経度BL（旧日本測地系）に換算しテキストボックスに書き込む
CalcXYToBL(0, inLat, inLon, inX, inY); // Bessel楕円体=日本測地系
// テキストボックスのBL（旧日本測地系）を座標変換しBL（JGD2000系）としテキストボックスみ書き込む
ForwardInterpol(convflag, inLat, inLon);
// 「『地域毎のﾊﾟﾗﾒｰﾀ』で変換する。そのパラメータがない所は変換しない｡」を選択し，そのパラメータがなかったので変換されなかった場合は，これ以上計算しない。変換結果のテキストボックスに既に”-9999.”が入っている。
if (convflag == true)
// 緯度経度BL（JGD2000系）をテキストボックスから読み込み平面直角座標値XY（JGD2000系）に換算しテキストボックスに書き込む
CalcBLToXY(1, inX, inY);// GRS80楕円体=世界測地系
}
/// <summary>
/// XYからXYへの逆方向（日本測地系2000から旧日本測地系へ）の変換を行います[doXyInterpol2]
/// </summary>
/// <param name="convflag">地域ごとのパラメータで変換するか否かのフラグ</param>
/// <param name="inLat">入力緯度</param>
/// <param name="inLon">入力経度</param>
/// <param name="inX">入力x</param>
/// <param name="inY">入力y</param>
public void ReverseXYInterpol(bool convflag, double inLat, double inLon, double inX, double inY)
{
// 平面直角座標値XY（JGD2000系）をテキストボックスから読み込み緯度経度BL（JGD2000系）に換算しテキストボックスに書き込む
CalcXYToBL(2, inLat, inLon, inX, inY); // GRS80楕円体=世界測地系
// テキストボックスのBL（JGD2000系）を座標変換しBL（旧日本測地系）としテキストボックスみ書き込む
ReverseInterpol(convflag, inLat, inLon);
// 「『地域毎のﾊﾟﾗﾒｰﾀ』で変換する。そのパラメータがない所は変換しない｡」を選択し，そのパラメータがなかったので変換されなかった場合は，これ以上計算しない。
// 変換結果のテキストボックスに既に”-9999.”が入っている。
if (convflag == true)
// 緯度経度BL（旧日本測地系）をテキストボックスから読み込み平面直角座標値XY（旧日本測地系）に換算しテキストボックスに書き込む
CalcBLToXY(1, inX, inY);// Bessel楕円体=日本測地系
}
/// <summary>
/// 「精密測地網一次基準点測量計算式」を導入し多くの項を採用。[doCalcBl2xyFile]
/// TKY2JGDではGRS80楕円体専門で利用。とはいってもどちらでも使える。called by XYFile
/// doCalcBl2xyのﾌｧｲﾙ一括処理用の値渡し版。ただし，ファイルのオープン，読み，書き，クローズは行わない。
/// </summary>
/// <param name="preCalc">True（最初に呼ばれるとき）なら関係変数を計算しstatic変数に代入し，Falseなら計算しない。</param>
/// <param name="ellipSuffix">楕円体配列の添え字：0=Bessel楕円体=日本測地系, 1=GRS80楕円体=世界測地系</param>
/// <param name="inB">緯度 単位は度</param>
/// <param name="inL">経度 単位は度</param>
/// <param name="outX">出力する平面直角座標X</param>
/// <param name="outY">出力する平面直角座標Y</param>
/// <param name="inLat">緯度</param>
/// <param name="inLon">経度</param>
public void CalcBLToXYFile(bool preCalc, int ellipSuffix, double inB, double inL, double outX, double outY, double inLat, double inLon)
{
double meridianCoefficiant = 0.0; //系番号，基準子午線の縮尺係数
double originBeta = 0.0; //原点の緯度，経度。基本的にradian
double originLambda = 0.0; //原点の緯度，経度。基本的にradian
double inBetaRadian = 0.0; // 入力される緯度，経度。基本的にradian
double inLambdaRadian = 0.0; // 入力される緯度，経度。基本的にradian
int deg = 0;
int menuite = 0;
double second = 0.0;
int sign = 0;
double gamma = 0.0; // =γ 子午線収差角。radian
double meridianScale = 0.0; // 縮尺係数
if (preCalc)
{
// 基準子午線の縮尺係数
meridianCoefficiant = 0.9999;
// コンボボックスから座標系，座標系原点緯度，経度の読み取り
double cB1 = inLat; //原点の緯度，経度。c = combo box
double cL1 = inLon; //原点の緯度，経度。c = combo box
DmsToDeg(cB1, ref originBeta, ref sign, ref deg, ref menuite, ref second);
DmsToDeg(cL1, ref originLambda, ref sign, ref deg, ref menuite, ref second);
originBeta = originBeta * DegToRad;
originLambda = originLambda * DegToRad; // 座標系原点の経度(rad)
}
try
{
inBetaRadian = inB * DegToRad;
inLambdaRadian = inL * DegToRad; // 変換したい経度(rad)
// 本計算 EP(0)はBessel楕円体=日本測地系，EP(1)はGRS-80楕円体=世界測地系
CalcBLToXY(preCalc, inBetaRadian, inLambdaRadian, originBeta, originLambda, meridianCoefficiant, ellipObj[ellipSuffix], ref outX, ref outY, ref gamma, ref meridianScale);
return; // エラー処理ルーチンが実行されないように Sub を終了します。
}
catch (OverflowException )
{
CommonUtils.LogWrite(string.Format(ConstString.Error.FailedLatLon, MethodBase.GetCurrentMethod().Name, inBetaRadian + ":" + inLambdaRadian));
return;
}
catch (Exception e)
{
CommonUtils.LogWrite(string.Format(ConstString.Error.ExceptionError, MethodBase.GetCurrentMethod().Name, e.Message));
return;
}
}
/// <summary>
/// 緯度,経度を平面直角座標X,Yに換算します。[doCalcBl2xy]
/// 「精密測地網一次基準点測量計算式」を導入し多くの項を採用。
/// </summary>
/// <param name="ellipSuffix">楕円体配列の添え字：0=Bessel楕円体=日本測地系, 1=GRS80楕円体=世界測地系</param>
/// <param name="inX">平面直角座標値，原点での北向き成分(meter)</param>
/// <param name="inY">平面直角座標値，原点での東向き成分(meter)</param>
public void CalcBLToXY(int ellipSuffix, double inX, double inY)
{
double meridianCoefficiant = 0.0; // 系番号，基準子午線の縮尺係数 ex.0.9999(X,Y), 0.9996(UTM)
double originBeta = 0.0; // 原点の緯度，経度。基本的にradian
double originLambda = 0.0; // 原点の緯度，経度。基本的にradian
double inBetaRadian = 0.0; // 入力される緯度，経度。基本的にradian
double inLambdaRadian = 0.0; // 入力される緯度，経度。基本的にradian
int degree = 0;
int minute = 0;
double second = 0.0;
int sign = 0;
string signStr = "";
double x = 0.0; // 求められる平面直角座標X,Y
double y = 0.0; // 求められる平面直角座標X,Y
double gamma = 0.0; // =γ 子午線収差角。radian
double gammaDeg = 0.0; // 真北方向角。deg
double meridianScale = 0.0; // 縮尺係数
// この計算は変換の第３段階なので緯度経度入力値のﾁｪｯｸは不要
// 基準子午線の縮尺係数
meridianCoefficiant = 0.9999;
try
{
// 入力
if (ellipSuffix == 1)
{
DmsToDeg(inX, ref inBetaRadian, ref sign, ref degree, ref minute, ref second);
DmsToDeg(inY, ref inLambdaRadian, ref sign, ref degree, ref minute, ref second);
}
else
{
DmsToDeg(inX, ref inBetaRadian, ref sign, ref degree, ref minute, ref second);
DmsToDeg(inY, ref inLambdaRadian, ref sign, ref degree, ref minute, ref second);
}
// degからradianへ換算
inBetaRadian = inBetaRadian * DegToRad; // 変換したい緯度(rad)
inLambdaRadian = inLambdaRadian * DegToRad; // 変換したい経度(rad)
// コンボボックスから座標系，座標系原点緯度，経度(in DMS format)の読み取り
double cB1 = inX;
double cL1 = inY;
DmsToDeg(cB1, ref originBeta, ref sign, ref degree, ref minute, ref second);
DmsToDeg(cL1, ref originLambda, ref sign, ref degree, ref minute, ref second);
originBeta = originBeta * DegToRad; // 座標系原点の緯度(rad)
originLambda = originLambda * DegToRad; // 座標系原点の経度(rad)
// 本計算 EP(0)はBessel楕円体=日本測地系，EP(1)はGRS-80楕円体=世界測地系
// True：楕円体，座標系に関わる変数の計算を行う
CalcBLToXY(true, inBetaRadian, inLambdaRadian, originBeta, originLambda, meridianCoefficiant, ellipObj[ellipSuffix], ref x, ref y, ref gamma, ref meridianScale);
// Ｘ，Ｙの計算 in meter
string outX = $"{x:######.0000}";
string outY = $"{y:######.0000}";
// 子午線収差角Gamma(rad)から，真北方向角(degree)へ
gammaDeg = -(gamma * RadToDeg);
DegToDmsNotZero3(gammaDeg, ref degree, ref minute, ref second, ref signStr);
string outDeg = $"{signStr}{degree:###}{minute:00}{second:00.000}";
// 縮尺係数 cf.旧成果表は小数点以下６桁
string outParam = $"{meridianScale:#.00000000}";
// outX,outY,outDeg,outParam に出力データが格納されているので、呼び出しもとで使えるように調整を。
return;
}
catch( OverflowException )
{
CommonUtils.LogWrite(string.Format(ConstString.Error.FailedLatLon, MethodBase.GetCurrentMethod().Name, inBetaRadian + ":" + inLambdaRadian));
return;
}
catch (Exception e)
{
CommonUtils.LogWrite(string.Format(ConstString.Error.ExceptionError, MethodBase.GetCurrentMethod().Name, e.Message));
return;
}
}
/// <summary>
/// called by picCalc1_click, calls Bilinear[doInterpol1]
/// 順方向（日本測地系からJGD2000）変換。"1"は順方向の意味。
/// （１）テキストボックスからデータを読み込み，
/// （２）パラメータを参照し，bilinear補間法により補間し，
/// （３）補間した，dL,dBを実用成果のL,Bに加える。
/// （４）結果をテキストボックスに書き込む。
/// </summary>
/// <param name="convflag">地域ごとのパラメータで変換するか否かのフラグ</param>
/// <param name="inBeta">緯度</param>
/// <param name="inLambda">経度</param>
public void ForwardInterpol(bool convflag, double inBeta, double inLambda)
{
int sign = 0;
string signStr = "";
int degree = 0;
int minute = 0;
double second = 0.0;
// Y:緯度，X:経度 （度）
double x1 = 0.0;
double y1 = 0.0;
double x2 = 0.0;
double y2 = 0.0;
double diffLambda = 0.0;
double diffBeta = 0.0;
TranPara tranPara = new TranPara();
//出力系
string outStatus = "";
string outLat = "";
string outLon = "";
// （１）テキストボックスからデータを読み込み，度分秒形式を度にし，
DmsToDeg(inBeta, ref y1, ref sign, ref degree, ref minute, ref second);
// 日本の領土外ならMessageBoxを表示
// 緯度のチェック
if (sign < 0 || degree < 20 || degree > 46)
MessageBox.Show(
string.Format(ConstString.Warning.LatLonOutofJapan, $"{(sign * degree):#0度}{minute:00分}{second:00.###秒}"),
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Warning
);
DmsToDeg(inLambda, ref x1, ref sign, ref degree, ref minute, ref second);
//経度のチェック
if (sign < 0 || degree < 120 || degree > 154)
MessageBox.Show(
string.Format(ConstString.Warning.LatLonOutofJapan, $"{(sign * degree):#0度}{minute:00分}{second:00.###秒}"),
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Warning
);
// (A) もし，「すべて3ﾊﾟﾗﾒｰﾀで変換」なら
if (allThreeParamFlag)
{
Tokyo97toITRF94(y1, x1, ref y2, ref x2);
outStatus = "「東京大正」三角点における「3ﾊﾟﾗﾒｰﾀ」で変換しました。\r\n日本測地系から世界測地系へ（→右方向に）変換しました。";
}
else
{
//「地域毎のﾊﾟﾗﾒｰﾀで変換」を試みる
// （２）パラメータを参照し，bilinear補間法により補間し，
tranPara = Bilinear1(y1, x1);
diffBeta = tranPara.DiffBeta;
diffLambda = tranPara.DiffLambda;
if (diffLambda == -9999.0 || diffBeta == -9999.0)
{
//「地域毎のﾊﾟﾗﾒｰﾀ」がなかったら
if (inputOptR3Flag)
{
//「3ﾊﾟﾗﾒｰﾀ」で変換する
Tokyo97toITRF94(y1, x1, ref y2, ref x2);
convflag = true;
outStatus = "「地域毎の変換ﾊﾟﾗﾒｰﾀ」がなかったので，「東京大正」三角点における「3ﾊﾟﾗﾒｰﾀ」で変換しました。\r\n日本測地系から世界測地系へ（→右方向に）変換しました。";
}
else
// 変換はしない。「『地域毎のﾊﾟﾗﾒｰﾀ』で変換する。そのパラメータがない所は変換しない｡」を選択した場合の処理。
{
// 変換はしないログ出力して終了
CommonUtils.LogWrite(ConstString.Warning.NoConvert);
return;
}
}
else
{
//「地域毎のﾊﾟﾗﾒｰﾀ」があったので変換する
// （３）補間した，dL,dBを実用成果のL,Bに加える
x2 = x1 + diffLambda / 3600.0;
y2 = y1 + diffBeta / 3600.0;
convflag = true;
outStatus = "「地域毎の変換ﾊﾟﾗﾒｰﾀ」で変換しました。\r\n日本測地系から世界測地系へ（→右方向に）変換しました。";
}
}
// （４）結果をテキストボックスに書き込む。
DegToDms(y2, ref degree, ref minute, ref second, ref signStr);
outLat = $"{signStr}{degree:#0}{minute:00}{second:00.00000}";
DegToDms(x2, ref degree, ref minute, ref second, ref signStr);
outLon = $"{signStr}{degree:#0}{minute:00}{second:00.00000}";
// outStaus, outLat, outLon に出力値が入っているので呼び出し下で受け取れるように調整を。
return;
}
/// <summary>
/// called by picCalc2_click, calls Bilinear
/// 逆方向（JGD2000から日本測地系）変換。"2"は逆方向の意味。
/// ２段階によるIteration(反復計算)により，パラメータファイル参照の緯度経度を日本測地系にする
/// （１）テキストボックスからデータを読み込み，
/// （２）パラメータを参照し，bilinear補間法により補間し，
/// （３）補間した，dL,dBを実用成果のL,Bに加える。
/// （４）結果をテキストボックスに書き込む。
/// </summary>
/// <param name="convflag">地域ごとのパラメータで変換するか否かのフラグ</param>
/// <param name="inBeta">緯度</param>
/// <param name="inLabmda">経度</param>
public void ReverseInterpol(bool convflag, double inBeta, double inLabmda)
{
int sign = 0;
string signStr = "";
int degree = 0;
int minute = 0;
double second = 0.0;
double x1 = 0.0;
double y1 = 0.0;
double x2 = 0.0;
double y2 = 0.0;
double diffLambda = 0.0;
double diffBeta = 0.0;
double tempLambda = 0.0;
double tempBeta = 0.0;
//出力系
string outStatus = "";
string outLat = "";
string outLon = "";
// （１）テキストボックスからデータを読み込み，度分秒形式を度にし，
DmsToDeg(inBeta, ref y1, ref sign, ref degree, ref minute, ref second);
// 日本の領土外ならMessageBoxを表示
// 緯度のチェック
if (sign < 0 || degree < 20 || degree > 46)
MessageBox.Show(
string.Format(ConstString.Warning.LatLonOutofJapan, $"{(sign * degree):#0度}{minute:00分}{second:00.###秒}"),
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Warning
);
DmsToDeg(inLabmda, ref x1, ref sign, ref degree, ref minute, ref second);
if (sign < 0 || degree < 120 || degree > 154)
MessageBox.Show(
string.Format(ConstString.Warning.LatLonOutofJapan, $"{(sign * degree):#0度}{minute:00分}{second:00.###秒}"),
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Warning
);
// (A) もし，「すべて3ﾊﾟﾗﾒｰﾀで変換」なら
if (allThreeParamFlag)
{
ITRF94toTokyo97(y1, x1, ref y2, ref x2);
convflag = true;
outStatus = "「東京大正」三角点における「3ﾊﾟﾗﾒｰﾀ」で変換しました。\r\n世界測地系から日本測地系へ（←左方向に）変換しました。";
}
else
{
//「地域毎のﾊﾟﾗﾒｰﾀで変換」を試みる
// （２）パラメータを参照し，bilinear補間法により補間し，
// (2-1)第一段階｛パラメータを参照しにいく緯度，経度は仮の値である。｝
tempLambda = x1 + 12.0 / 3600.0; // 平均値を用いてJGD2000系を日本測地系へ
tempBeta = y1 - 12.0 / 3600.0; // 平均値を用いてJGD2000系を日本測地系へ
TranPara tranPara = Bilinear1(tempBeta, tempLambda);
diffBeta = tranPara.DiffBeta;
diffLambda = tranPara.DiffLambda;
if (diffLambda == -9999.0 || diffBeta == -9999.0)
{
//「地域毎のﾊﾟﾗﾒｰﾀ」がなかったら
if (inputOptR3Flag)
{
//「3ﾊﾟﾗﾒｰﾀ」で変換する
ITRF94toTokyo97(y1, x1, ref y2, ref x2);
convflag = true;
outStatus = "「地域毎の変換ﾊﾟﾗﾒｰﾀ」がなかったので，「東京大正」三角点における「3ﾊﾟﾗﾒｰﾀ」で変換しました。\r\n世界測地系から日本測地系へ（←左方向に）変換しました。";
}
else
{
// 変換はしないログ出力して終了
CommonUtils.LogWrite(ConstString.Warning.NoConvert);
return;
}
}
else
{
//「地域毎のﾊﾟﾗﾒｰﾀ」があったので変換する
// （３）補間した，-dL,-dBを実用成果のL,Bに加える
// (3-1)第一段階
tempLambda = x1 - diffLambda / 3600.0;
tempBeta = y1 - diffBeta / 3600.0;
// (2-2)第二段階｛パラメータを参照しにいく緯度，経度はほぼ真の日本測地系の値である。｝
tranPara = Bilinear1(tempBeta, tempLambda);
diffBeta = tranPara.DiffBeta;
diffLambda = tranPara.DiffLambda;
if (diffLambda == -9999.0 || diffBeta == -9999.0)
{
// 変換はしないログ出力して終了
CommonUtils.LogWrite(ConstString.Warning.NoConvert);
return;
}
// （３）補間した，-dL,-dBを実用成果のL,Bに加える
// (3-1)第二段階
x2 = x1 - diffLambda / 3600.0;
y2 = y1 - diffBeta / 3600.0;
convflag = true;
outStatus = "「地域毎の変換ﾊﾟﾗﾒｰﾀ」で変換しました。\r\n世界測地系から日本測地系へ（←左方向に）変換しました。";
} // 「地域毎のﾊﾟﾗﾒｰﾀ」があるかないかの判断おしまい
} // 「3ﾊﾟﾗﾒｰﾀ」か「地域毎のﾊﾟﾗﾒｰﾀ」かの判断おしまい
// （４）結果をテキストボックスに書き込む。
DegToDms(y2, ref degree, ref minute, ref second, ref signStr);
outLat = $"{signStr}{degree:#0}{minute:00}{second:00.00000}";
DegToDms(x2, ref degree, ref minute, ref second, ref signStr);
outLon = $"{signStr}{degree:#0}{minute:00}{second:00.00000}";
// outStaus, outLat, outLon に出力値が入っているので呼び出し下で受け取れるように調整を。
return;
}
/// <summary>
/// DMS(DDDMMSS.SSS)をdegreeにする。
/// </summary>
/// <param name="inDMS">DMS形式の入力座標</param>
/// <param name="deg">度数表示の出力</param>
/// <param name="sign">符号 + or -</param>
/// <param name="degree">度</param>
/// <param name="minute">分</param>
/// <param name="second">秒</param>
public void DmsToDeg(double inDMS, ref double deg, ref int sign, ref int degree, ref int minute, ref double second)
{
sign = Math.Sign(inDMS);
// N88BASICの場合,小さな数を足さないと、1000.つまり10'を入力したとき、9'99.99"
// が入力されたと解釈されてしまい、40秒のずれが生じてしまう。
// Quick BASICとVisual Basicでは大丈夫。以上，1991,1992年
// しかし、1992/9/29の調査によると，QBでも、90100.を入力すると、分の計算で，
// 0.010と解釈してほしいところ，0.00999999999...と解釈されることから，
// 以下のように、小さい数を加える。このこれより大きくても小さくてもなかなか
// うまくいかい。この問題は，100.の場合は起こらない。 飛田幹男 >>Ver.2.1
// また，7/18/1993 TRN93作成時に、334000. を入力するとだめ。そこで、Ver.2.2.
double tmpDMS = Math.Abs(inDMS / 10000.0) + 0.0000000000001;
degree = (int)tmpDMS;
double absoluteMS = (tmpDMS - degree) * 100.0;
minute = (int)absoluteMS;
second = (absoluteMS - minute) * 100.0;
deg = (degree + (minute + second / 60.0) / 60.0) * sign;
}
/// <summary>
/// degreeをDMS(DDDMMSS.SSS)にする。
/// このサブプログラムは，秒の小数点以下「５桁」の場合に「60.00000」秒にならないようになっている。
/// deg(度)を符号SNSとD,AM,Sにする。
/// </summary>
/// <param name="inDegree">緯度、経度(arcsec)</param>
/// <param name="degree">度</param>
/// <param name="minute">分</param>
/// <param name="second">秒</param>
/// <param name="signStr">符号文字列 " " or "-"</param>
public void DegToDms(double inDegree, ref int degree, ref int minute, ref double second, ref string signStr)
{
if ( Math.Sign(inDegree) < 0)
signStr = "-";
else
signStr = " ";
double adeg = Math.Abs(inDegree) + 0.00000000001;
degree = (int)adeg;
double absoluteMinute = (adeg - degree) * 60.0;
minute = (int)absoluteMinute;
second = (absoluteMinute - minute) * 60.0;
// 注意：下のFormat括弧内の"00.000000"の小数点以下の0の数は実際に出力するときの数に合わせること。
// 1495960.000000を避け，1500000.000000になるように
if ( $"{second:00.00000}".Substring(0, 2) == "60" )
{
second = 0.0;
minute = minute + 1;
if (minute == 60)
{
minute = 0;
degree = degree + 1;
}
}
}
/// <summary>
/// degreeをDMS(DDDMMSS.SSS)にする。
/// このサブプログラムは，秒の小数点以下「３桁」の場合に「60.000」秒にならないようになっている。
/// deg(度)を符号SNSとD,AM,Sにする。
/// </summary>
/// <param name="inDeg">緯度、経度(arcsec)</param>
/// <param name="degree">度</param>
/// <param name="minute">分</param>
/// <param name="second">秒</param>
/// <param name="signStr">符号 " " or "-"</param>
public void DegToDmsNotZero3(double inDeg, ref int degree, ref int minute, ref double second, ref string signStr)
{
int sign = Math.Sign(inDeg);
if (sign < 0)
signStr = "-";
else
signStr = " ";
double adeg = Math.Abs(inDeg) + 0.00000000001;
degree = (int)adeg;
double absoluteMinute = (adeg - degree) * 60.0;
minute = (int)absoluteMinute;
second = (absoluteMinute - minute) * 60.0;
// 注意：下のFormat括弧内の"00.000000"の小数点以下の0の数は実際に出力するときの数に合わせること。
// 1495960.000000を避け，1500000.000000になるように
if ( $"{second:00.000}".Substring(0, 2) == "60")
{
second = 0.0;
minute = minute + 1;
if (minute == 60)
{
minute = 0;
degree = degree + 1;
}
}
}
/// <summary>
/// 緯度,経度を平面直角座標X,Yに換算します。called by doCalcBl2xy,doCalcBl2xyFile
/// 「精密測地網一次基準点測量計算式」を導入し多くの項を採用。
/// </summary>
/// <param name="preCalc"></param>
/// <param name="inBeta">入力された緯度</param>
/// <param name="inLambda">入力された経度</param>
/// <param name="beta">座標系原点の緯度</param>
/// <param name="lambda">座標系原点の経度</param>
/// <param name="meridianCoefficiant">基準子午線の縮尺係数 ex.0.9999(X,Y), 0.9996(UTM)</param>
/// <param name="EP">楕円体のパラメータ入り構造体</param>
/// <param name="X">平面直角座標値，原点での北向き成分(meter)</param>
/// <param name="Y">平面直角座標値，原点での東向き成分(meter)</param>
/// <param name="gamma">子午線収差角(radian),これのマイナスが真北方向角</param>
/// <param name="outMeridian">縮尺係数の出力</param>
private void CalcBLToXY(
bool preCalc,
double inBeta,
double inLambda,
double beta,
double lambda,
double meridianCoefficiant,
EllipParam EP,
ref double X,
ref double Y,
ref double gamma,
ref double outMeridian)
{
double AEE = 0.0;
double CEE = 0.0;
double Ep2 = 0.0;
double e2 = 0.0;
double e4 = 0.0;
double e6 = 0.0;
double e8 = 0.0;
double e10 = 0.0;
double e12 = 0.0;
double e14 = 0.0;
double e16 = 0.0;
double AJ = 0.0;
double BJ = 0.0;
double CJ = 0.0;
double DJ = 0.0;
double EJ = 0.0;
double FJ = 0.0;
double GJ = 0.0;
double HJ = 0.0;
double IJ = 0.0;
double S0 = 0.0;
double dL = 0.0;
double DL2 = 0.0;
double DL4 = 0.0;
double DL6 = 0.0;
double S = 0.0; // 赤道から求点までの子午線長の計算
double COS2 = 0.0;
double Eta2phi = 0.0;
double Mphi = 0.0;
double Nphi = 0.0; // phi(=B)の関数
double T = 0.0;
double T2 = 0.0;
double T4 = 0.0;
double T6 = 0.0;
dL = inLambda - lambda; // Δλ
if (preCalc)
{
e2 = EP.Eccentricity;
e4 = e2 * e2;
e6 = e4 * e2;
e8 = e4 * e4;
e10 = e8 * e2;
e12 = e8 * e4;
e14 = e8 * e6;
e16 = e8 * e8;
// 定数項
AEE = EP.Axix * (1.0 - EP.Eccentricity); // a(1-e2)
CEE = EP.Axix / (double)Math.Sqrt(1.0 - EP.Eccentricity); // C=a*Math.Sqrt(1+e'2)=a/Math.Sqrt(1-e2)
Ep2 = EP.Eccentricity / (double)(1.0 - EP.Eccentricity); // e'2 (e prime 2) Eta2phiを計算するため
// 「緯度を与えて赤道からの子午線弧長を求める計算」のための９つの係数を求める。
// 「精密測地網一次基準点測量計算式」P55,P56より。係数チェック済み1999/10/19。
AJ = 4927697775.0 / 7516192768.0 * e16;
AJ = AJ + 19324305.0 / 29360128.0 * e14;
AJ = AJ + 693693.0 / 1048576.0 * e12;
AJ = AJ + 43659.0 / 65536.0 * e10;
AJ = AJ + 11025.0 / 16384.0 * e8;
AJ = AJ + 175.0 / 256.0 * e6;
AJ = AJ + 45.0 / 64.0 * e4;
AJ = AJ + 3.0 / 4.0 * e2;
AJ = AJ + 1.0;
BJ = 547521975.0 / 469762048.0 * e16;
BJ = BJ + 135270135.0 / 117440512.0 * e14;
BJ = BJ + 297297.0 / 262144.0 * e12;
BJ = BJ + 72765.0 / 65536.0 * e10;
BJ = BJ + 2205.0 / 2048.0 * e8;
BJ = BJ + 525.0 / 512.0 * e6;
BJ = BJ + 15.0 / 16.0 * e4;
BJ = BJ + 3.0 / 4.0 * e2;
CJ = 766530765.0 / 939524096.0 * e16;
// CJ = CJ + 45090045.0 / 5870256.0 * e14 精密測地網一次基準点測量作業規定の誤りによるバグ
CJ = CJ + 45090045.0 / 58720256.0 * e14;
CJ = CJ + 1486485.0 / 2097152.0 * e12;
CJ = CJ + 10395.0 / 16384.0 * e10;
CJ = CJ + 2205.0 / 4096.0 * e8;
CJ = CJ + 105.0 / 256.0 * e6;
CJ = CJ + 15.0 / 64.0 * e4;
DJ = 209053845.0 / 469762048.0 * e16;
DJ = DJ + 45090045.0 / 117440512.0 * e14;
DJ = DJ + 165165.0 / 524288.0 * e12;
DJ = DJ + 31185.0 / 131072.0 * e10;
DJ = DJ + 315.0 / 2048.0 * e8;
DJ = DJ + 35.0 / 512.0 * e6;
EJ = 348423075.0 / 1879048192.0 * e16;
EJ = EJ + 4099095.0 / 29360128.0 * e14;
EJ = EJ + 99099.0 / 1048576.0 * e12;
EJ = EJ + 3465.0 / 65536.0 * e10;
EJ = EJ + 315.0 / 16384.0 * e8;
FJ = 26801775.0 / 469762048.0 * e16;
FJ = FJ + 4099095.0 / 117440512.0 * e14;
FJ = FJ + 9009.0 / 524288.0 * e12;
FJ = FJ + 693.0 / 131072.0 * e10;
GJ = 11486475.0 / 939524096.0 * e16;
GJ = GJ + 315315.0 / 58720256.0 * e14;
GJ = GJ + 3003.0 / 2097152.0 * e12;
HJ = 765765.0 / 469762048.0 * e16;
HJ = HJ + 45045.0 / 117440512.0 * e14;
IJ = 765765.0 / 7516192768.0 * e16;
// 赤道から座標系原点までの子午線長の計算
MeridS(beta, AEE, AJ, BJ, CJ, DJ, EJ, FJ, GJ, HJ, IJ, ref S0);
}
// 何度も使う式を変数に代入
T = Math.Tan(inBeta);
T2 = T * T; T4 = T2 * T2; T6 = T4 * T2;
COS2 = Math.Cos(inBeta) * Math.Cos(inBeta);
Eta2phi = Ep2 * COS2; // =η1*η1
Mphi = CEE / Math.Sqrt(Math.Pow((1.0 + Eta2phi), 3.0)); // 「精密測地網一次基準点測量計算式」P52のM(phi)
Nphi = CEE / Math.Sqrt(1.0 + Eta2phi); // 「精密測地網一次基準点測量計算式」P52のN(phi)
DL2 = dL * dL; DL4 = DL2 * DL2; DL6 = DL4 * DL2;
// 赤道から求点までの子午線長の計算
MeridS(inBeta, AEE, AJ, BJ, CJ, DJ, EJ, FJ, GJ, HJ, IJ, ref S);
// Ｘ，Ｙの計算 in meter
// 「精密測地網一次基準点測量計算式」P52,53のx,yを求める式より
X = -(-1385.0 + 3111.0 * T2 - 543.0 * T4 + T6) * DL6 * Math.Pow(COS2, 3.0) / 40320.0;
X = X - (-61.0 + 58.0 * T2 - T4 - 270.0 * Eta2phi + 330.0 * T2 * Eta2phi) * DL4 * Math.Pow(COS2, 2.0) / 720.0;
X = X + (5.0 - T2 + 9.0 * Eta2phi + 4.0 * Eta2phi * Eta2phi) * DL2 * COS2 / 24.0;
X = X + 1.0 / 2.0;
X = X * Nphi * COS2 * T * DL2;
X = X + S - S0;
X = X * meridianCoefficiant;
Y = -(-61.0 + 479.0 * T2 - 179.0 * T4 + T6) * DL6 * Math.Pow(COS2, 3.0) / 5040.0;
Y = Y - (-5.0 + 18.0 * T2 - T4 - 14.0 * Eta2phi + 58.0 * T2 * Eta2phi) * DL4 * Math.Pow(COS2, 2.0) / 120.0;
Y = Y - (-1.0 + T2 - Eta2phi) * DL2 * COS2 / 6.0;
Y = Y + 1.0;
Y = Y * Nphi * Math.Cos(inBeta) * dL;
Y = Y * meridianCoefficiant;
// 子午線収差角 これのマイナスが真北方向角(rad)
// 「精密測地網一次基準点測量計算式」P53のγを求める式より
gamma = COS2 * COS2 * (2.0 - T2) * Math.Pow(dL, 4.0) / 15.0;
gamma = gamma + COS2 * (1.0 + 3.0 * Eta2phi + 2.0 * Math.Pow(Eta2phi, 2.0)) * Math.Pow(dL, 2.0) / 3.0;
gamma = gamma + 1.0;
gamma = gamma * Math.Cos(inBeta) * T * dL;
// 縮尺係数の計算 「精密測地網一次基準点測量計算式」P51のmを求める式より
outMeridian = Math.Pow(Y, 4.0) / (24.0 * Mphi * Mphi * Nphi * Nphi * Math.Pow(meridianCoefficiant, 4.0));
outMeridian = outMeridian + Y * Y / (2.0 * Mphi * Nphi * Math.Pow(meridianCoefficiant, 2.0));
outMeridian = outMeridian + 1.0;
outMeridian = outMeridian * meridianCoefficiant;
}
/// <summary>
/// X,YをB,Lに換算します。called by XYFile
/// このサブプログラムが複数回呼ばれることを想定すると，楕円体パラメータと座標系が変更されずに
/// 呼ばれた場合，2回目以降，関係変数をダブって計算するのは無駄なので，
/// PreCalcがTrue（最初に呼ばれるとき）なら関係変数を計算しstatic変数に代入し，Falseなら計算しない。
/// </summary>
/// <param name="PreCalc">既に計算済みの場合、false</param>
/// <param name="ellipSuffix">楕円体配列の添え字：0=Bessel楕円体=日本測地系, 1=GRS80楕円体=世界測地系</param>
/// <param name="X">平面直角座標</param>
/// <param name="Y">平面直角座標</param>
/// <param name="Bout">度単位の緯度</param>
/// <param name="Lout">度単位の経度</param>
/// <param name="inX">平面直角座標 入力値</param>
/// <param name="inY">平面直角座標 入力値</param>
private void DoCalcXy2blFile(bool PreCalc, int ellipSuffix, double X, double Y, ref double Bout, ref double Lout, double inX, double inY)
{
const double M0 = 0.9999; //系番号，基準子午線の縮尺係数
double cB1 = 0.0; //原点の緯度，経度。c = combo box
double cL1 = 0.0; //原点の緯度，経度。c = combo box
double B1=0.0; //原点の緯度，経度。基本的にradian
double L1 = 0.0; //原点の緯度，経度。基本的にradian
double AEE = 0.0;
double CEE = 0.0;
double Ep2 = 0.0;
EllipParam EPs = new EllipParam();
double AJ = 0.0;
double BJ = 0.0;
double CJ = 0.0;
double DJ = 0.0;
double EJ = 0.0;
double FJ = 0.0;
double GJ = 0.0;
double HJ = 0.0;
double IJ = 0.0;
double e2 = 0.0;
double e4 = 0.0;
double e6 = 0.0;
double e8 = 0.0;
double e10 = 0.0;
double e12 = 0.0;
double e14 = 0.0;
double e16 = 0.0;
double B = 0.0; // 求める緯度，経度。基本的にradian
double L = 0.0; // 求める緯度，経度。基本的にradian
int D = 0;
int AM = 0;
double S = 0.0;
int SN = 0;
double S0 = 0.0; // 赤道から座標原点までの子午線弧長
double M = 0.0;
double Eta2 = 0.0; // phi1の関数
double M1 = 0.0; // phi1の関数
double N1 = 0.0; // phi1の関数
double T = 0.0;
double T2 = 0.0;
double T4 = 0.0;
double T6 = 0.0;
double S1 = 0.0;
double phi1 = 0.0;
double oldphi1 = 0.0;
int icountphi1 = 0;
double Bunsi = 0.0;
double Bunbo = 0.0;
double YM0 = 0.0;
double N1CosPhi1 = 0.0;
if (PreCalc)
{
// 楕円体ﾊﾟﾗﾒｰﾀの代入
EPs.Axix = ellipObj[ellipSuffix].Axix;
EPs.Flattening1 = ellipObj[ellipSuffix].Flattening1;
EPs.Flattening = ellipObj[ellipSuffix].Flattening;
EPs.Eccentricity = ellipObj[ellipSuffix].Eccentricity;
EPs.Namec = ellipObj[ellipSuffix].Namec;
e2 = EPs.Eccentricity;
e4 = e2 * e2;
e6 = e4 * e2;
e8 = e4 * e4;
e10 = e8 * e2;
e12 = e8 * e4;
e14 = e8 * e6;
e16 = e8 * e8;
// 定数項 the same as bl2xy
AEE = EPs.Axix * (1.0 - EPs.Eccentricity); // a(1-e2)
CEE = EPs.Axix / (double)Math.Sqrt(1.0 - EPs.Eccentricity); // C=a*Math.Sqrt(1+e'2)=a/Math.Sqrt(1-e2)
Ep2 = EPs.Eccentricity / (double)(1.0 - EPs.Eccentricity); // e'2 (e prime 2) Eta2を計算するため
// 「緯度を与えて赤道からの子午線弧長を求める計算」のための９つの係数を求める。
// 「精密測地網一次基準点測量計算式」P55,P56より。係数チェック済み1999/10/19。
AJ = 4927697775.0 / 7516192768.0 * e16;
AJ = AJ + 19324305.0 / 29360128.0 * e14;
AJ = AJ + 693693.0 / 1048576.0 * e12;
AJ = AJ + 43659.0 / 65536.0 * e10;
AJ = AJ + 11025.0 / 16384.0 * e8;
AJ = AJ + 175.0 / 256.0 * e6;
AJ = AJ + 45.0 / 64.0 * e4;
AJ = AJ + 3.0 / 4.0 * e2;
AJ = AJ + 1.0;
BJ = 547521975.0 / 469762048.0 * e16;
BJ = BJ + 135270135.0 / 117440512.0 * e14;
BJ = BJ + 297297.0 / 262144.0 * e12;
BJ = BJ + 72765.0 / 65536.0 * e10;
BJ = BJ + 2205.0 / 2048.0 * e8;
BJ = BJ + 525.0 / 512.0 * e6;
BJ = BJ + 15.0 / 16.0 * e4;
BJ = BJ + 3.0 / 4.0 * e2;
CJ = 766530765.0 / 939524096.0 * e16;
// CJ = CJ + 45090045.0 / 5870256.0 * e14 精密測地網一次基準点測量作業規定の誤りによるバグ
CJ = CJ + 45090045.0 / 58720256.0 * e14;
CJ = CJ + 1486485.0 / 2097152.0 * e12;
CJ = CJ + 10395.0 / 16384.0 * e10;
CJ = CJ + 2205.0 / 4096.0 * e8;
CJ = CJ + 105.0 / 256.0 * e6;
CJ = CJ + 15.0 / 64.0 * e4;
DJ = 209053845.0 / 469762048.0 * e16;
DJ = DJ + 45090045.0 / 117440512.0 * e14;
DJ = DJ + 165165.0 / 524288.0 * e12;
DJ = DJ + 31185.0 / 131072.0 * e10;
DJ = DJ + 315.0 / 2048.0 * e8;
DJ = DJ + 35.0 / 512.0 * e6;
EJ = 348423075.0 / 1879048192.0 * e16;
EJ = EJ + 4099095.0 / 29360128.0 * e14;
EJ = EJ + 99099.0 / 1048576.0 * e12;
EJ = EJ + 3465.0 / 65536.0 * e10;
EJ = EJ + 315.0 / 16384.0 * e8;
FJ = 26801775.0 / 469762048.0 * e16;
FJ = FJ + 4099095.0 / 117440512.0 * e14;
FJ = FJ + 9009.0 / 524288.0 * e12;
FJ = FJ + 693.0 / 131072.0 * e10;
GJ = 11486475.0 / 939524096.0 * e16;
GJ = GJ + 315315.0 / 58720256.0 * e14;
GJ = GJ + 3003.0 / 2097152.0 * e12;
HJ = 765765.0 / 469762048.0 * e16;
HJ = HJ + 45045.0 / 117440512.0 * e14;
IJ = 765765.0 / 7516192768.0 * e16;
// 座標系番号，緯度，経度の読みとり
cB1 = inX; // 原点緯度
cL1 = inY; // 原点経度
DmsToDeg(cB1, ref B1, ref SN, ref D, ref AM, ref S);
DmsToDeg(cL1, ref L1, ref SN, ref D, ref AM, ref S);
B1 = B1 * DegToRad;
L1 = L1 * DegToRad;
}
try
{
// 赤道からの子午線長の計算
MeridS(B1, AEE, AJ, BJ, CJ, DJ, EJ, FJ, GJ, HJ, IJ, ref S0); // 赤道から座標原点までの子午線弧長
M = S0 + X / M0; // 「精密測地網一次基準点測量計算式」P57中段
// Baileyの式による異性緯度（isometric latitude）phi1の計算。
// 「精密測地網一次基準点測量計算式」P57の11.(1)の式から。
// この式と「現代測量学１ 測量の数学的基礎」P102の式とは，Math.Cos(phi1)だけ異なる。
// この式を導入したためベッセル楕円体以外で往復計算OKとなった。
icountphi1 = 0;
phi1 = B1;
do
{
icountphi1 = icountphi1 + 1;
oldphi1 = phi1;
MeridS(phi1, AEE, AJ, BJ, CJ, DJ, EJ, FJ, GJ, HJ, IJ, ref S1);
Bunsi = 2.0 * (S1 - M) * Math.Pow((1.0 - EPs.Eccentricity * Math.Sin(phi1) * Math.Sin(phi1)), 1.5);
Bunbo = 3.0 * EPs.Eccentricity * (S1 - M) * Math.Sin(phi1) * Math.Cos(phi1) * Math.Sqrt(1.0 - EPs.Eccentricity * Math.Sin(phi1) * Math.Sin(phi1)) - 2.0 * EPs.Axix * (1.0 - EPs.Eccentricity);
phi1 = phi1 + Bunsi / Bunbo;
}
while (!((Math.Abs(phi1 - oldphi1) < 0.00000000000001) || (icountphi1 > 100))); // 赤道から点までの子午線弧長
// 本では1e-12で十分 iterationの回数は４回
// 本計算
// 何度も使う式を変数に代入
YM0 = Y / M0;
T = Math.Tan(phi1); // 「精密測地網一次基準点測量計算式」P51のt1に等しい
T2 = T * T; T4 = T2 * T2; T6 = T4 * T2;
Eta2 = Ep2 * Math.Cos(phi1) * Math.Cos(phi1); // =η1*η1
M1 = CEE / Math.Sqrt(Math.Pow((1.0 + Eta2), 3.0));
N1 = CEE / Math.Sqrt(1.0 + Eta2);
N1CosPhi1 = N1 * Math.Cos(phi1);
// 緯度Bの計算 「精密測地網一次基準点測量計算式」P51のphiを求める式より
B = ((1385.0 + 3633.0 * T2 + 4095 * T4 + 1575.0 * T6) / (40320.0 * Math.Pow(N1, 8.0))) * Math.Pow(YM0, 8.0);
B = B - ((61.0 + 90.0 * T2 + 45 * T4 + 107.0 * Eta2 - 162.0 * T2 * Eta2 - 45.0 * T4 * Eta2) / (720.0 * Math.Pow(N1, 6.0))) * Math.Pow(YM0, 6.0);
B = B + ((5.0 + 3.0 * T2 + 6.0 * Eta2 - 6.0 * T2 * Eta2 - 3.0 * Math.Pow(Eta2, 2) - 9.0 * T2 * Math.Pow(Eta2, 2)) / (24.0 * Math.Pow(N1, 4.0))) * Math.Pow(YM0, 4.0);
B = B - ((1.0 + Eta2) / (2.0 * Math.Pow(N1, 2.0))) * Math.Pow(YM0, 2.0);
B = B * T;
B = B + phi1;
// 経度Lの計算 「精密測地網一次基準点測量計算式」P51のΔλを求める式より
L = -((61.0 + 662.0 * T2 + 1320.0 * T4 + 720.0 * T6) / (5040.0 * Math.Pow(N1, 6.0) * N1CosPhi1)) * Math.Pow(YM0, 7.0);
L = L + ((5.0 + 28.0 * T2 + 24.0 * T4 + 6.0 * Eta2 + 8.0 * T2 * Eta2) / (120.0 * Math.Pow(N1, 4.0) * N1CosPhi1)) * Math.Pow(YM0, 5.0);
L = L - ((1.0 + 2.0 * T2 + Eta2) / (6.0 * Math.Pow(N1, 2.0) * N1CosPhi1)) * Math.Pow(YM0, 3.0);
L = L + (1.0 / N1CosPhi1) * YM0;
L = L + L1;
// 縮尺係数と真北方向角の計算はファイル処理では不要
// 出力
Bout = B * RadToDeg;
Lout = L * RadToDeg;
return; // エラー処理ルーチンが実行されないように Sub を終了します。
}
catch(OverflowException oe)
{
CommonUtils.LogWrite(string.Format(ConstString.Error.FailedXY, MethodBase.GetCurrentMethod().Name, B,L));
return;
}
catch (Exception e)
{
CommonUtils.LogWrite(string.Format(ConstString.Error.FailedXY, MethodBase.GetCurrentMethod().Name, B, L));
return;
}
}
/// <summary>
/// 平面直角座標X,Yを緯度,経度に換算します。
/// 「精密測地網一次基準点測量計算式」を導入し多くの項を採用。
/// NAD83 Datum Sheet Dataのアラスカの点で往復計算の経度差が0.00156秒あったのを0.00000秒にした。
/// この計算は変換の第1段階なので緯度経度入力値のﾁｪｯｸが必要
/// </summary>
/// <param name="ellipSuffix">楕円体配列の添え字：0=Bessel楕円体=日本測地系, 1=GRS80楕円体=世界測地系</param>
/// <param name="inLat">入力緯度</param>
/// <param name="inLon">入力経度</param>
/// <param name="inX">入力X</param>
/// <param name="inY">入力Y</param>
private void CalcXYToBL(int ellipSuffix, double inLat, double inLon, double inX, double inY)
{
double X = 0.0; // 平面直角座標。入力値。meter
double Y = 0.0; // 平面直角座標。入力値。meter
double cB1 = 0.0; // 原点の緯度，経度。c = combo box
double cL1 = 0.0; // 原点の緯度，経度。c = combo box
double betaOrigin = 0.0; // 原点の緯度，経度。基本的にradian
double lambdaOrigin = 0.0; // 原点の緯度，経度。基本的にradian
double beta = 0.0; // 求める緯度，経度。基本的にradian
double lambda = 0.0; // 求める緯度，経度。基本的にradian
double gamma = 0.0; // =γ 子午線収差角。radian
int degree = 0;
int minute = 0;
double second = 0.0;
int sign = 0;
string signString = "";
double S0 = 0.0; // 赤道から座標原点までの子午線弧長
double M = 0.0;
double T4 = 0.0;
double T6 = 0.0;
double S1 = 0.0;
double phi1 = 0.0;
double oldphi1 = 0.0;
int icount = 0;
double Bunsi = 0.0;
double Bunbo = 0.0;
EllipParam EPs = new EllipParam();
//出力系
string outX = "";
string outY = "";
string outLat = "";
string outLon = "";
string outDeg = "";
string outParam = "";
// 基準子午線の縮尺係数
double M0 = 0.9999;
// 楕円体ﾊﾟﾗﾒｰﾀの代入
EPs.Axix = ellipObj[ellipSuffix].Axix;
EPs.Flattening1 = ellipObj[ellipSuffix].Flattening1;
EPs.Flattening = ellipObj[ellipSuffix].Flattening;
EPs.Eccentricity = ellipObj[ellipSuffix].Eccentricity;
EPs.Namec = ellipObj[ellipSuffix].Namec;
double e2 = EPs.Eccentricity;
double e4 = e2 * e2;
double e6 = e4 * e2;
double e8 = e4 * e4;
double e10 = e8 * e2;
double e12 = e8 * e4;
double e14 = e8 * e6;
double e16 = e8 * e8;
// 定数項 the same as bl2xy
double AEE = EPs.Axix * (1.0 - EPs.Eccentricity); // a(1-e2)
double CEE = EPs.Axix / (double)Math.Sqrt(1.0 - EPs.Eccentricity); // C=a*Math.Sqrt(1+e'2)=a/Math.Sqrt(1-e2)
double Ep2 = EPs.Eccentricity / (double)(1.0 - EPs.Eccentricity); // e'2 (e prime 2) Eta2を計算するため
// 「緯度を与えて赤道からの子午線弧長を求める計算」のための９つの係数を求める。
// 「精密測地網一次基準点測量計算式」P55,P56より。係数チェック済み1999/10/19。
double AJ = 4927697775.0 / 7516192768.0 * e16;
AJ = AJ + 19324305.0 / 29360128.0 * e14;
AJ = AJ + 693693.0 / 1048576.0 * e12;
AJ = AJ + 43659.0 / 65536.0 * e10;
AJ = AJ + 11025.0 / 16384.0 * e8;
AJ = AJ + 175.0 / 256.0 * e6;
AJ = AJ + 45.0 / 64.0 * e4;
AJ = AJ + 3.0 / 4.0 * e2;
AJ = AJ + 1.0;
double BJ = 547521975.0 / 469762048.0 * e16;
BJ = BJ + 135270135.0 / 117440512.0 * e14;
BJ = BJ + 297297.0 / 262144.0 * e12;
BJ = BJ + 72765.0 / 65536.0 * e10;
BJ = BJ + 2205.0 / 2048.0 * e8;
BJ = BJ + 525.0 / 512.0 * e6;
BJ = BJ + 15.0 / 16.0 * e4;
BJ = BJ + 3.0 / 4.0 * e2;
double CJ = 766530765.0 / 939524096.0 * e16;
// CJ = CJ + 45090045.0 / 5870256.0 * e14 精密測地網一次基準点測量作業規定の誤りによるバグ
CJ = CJ + 45090045.0 / 58720256.0 * e14;
CJ = CJ + 1486485.0 / 2097152.0 * e12;
CJ = CJ + 10395.0 / 16384.0 * e10;
CJ = CJ + 2205.0 / 4096.0 * e8;
CJ = CJ + 105.0 / 256.0 * e6;
CJ = CJ + 15.0 / 64.0 * e4;
double DJ = 209053845.0 / 469762048.0 * e16;
DJ = DJ + 45090045.0 / 117440512.0 * e14;
DJ = DJ + 165165.0 / 524288.0 * e12;
DJ = DJ + 31185.0 / 131072.0 * e10;
DJ = DJ + 315.0 / 2048.0 * e8;
DJ = DJ + 35.0 / 512.0 * e6;
double EJ = 348423075.0 / 1879048192.0 * e16;
EJ = EJ + 4099095.0 / 29360128.0 * e14;
EJ = EJ + 99099.0 / 1048576.0 * e12;
EJ = EJ + 3465.0 / 65536.0 * e10;
EJ = EJ + 315.0 / 16384.0 * e8;
double FJ = 26801775.0 / 469762048.0 * e16;
FJ = FJ + 4099095.0 / 117440512.0 * e14;
FJ = FJ + 9009.0 / 524288.0 * e12;
FJ = FJ + 693.0 / 131072.0 * e10;
double GJ = 11486475.0 / 939524096.0 * e16;
GJ = GJ + 315315.0 / 58720256.0 * e14;
GJ = GJ + 3003.0 / 2097152.0 * e12;
double HJ = 765765.0 / 469762048.0 * e16;
HJ = HJ + 45045.0 / 117440512.0 * e14;
double IJ = 765765.0 / 7516192768.0 * e16;
// 座標系番号，緯度，経度の読みとり
cB1 = inX; // 原点緯度
cL1 = inY; // 原点経度
DmsToDeg(cB1, ref betaOrigin, ref sign, ref degree, ref minute, ref second);
DmsToDeg(cL1, ref lambdaOrigin, ref sign, ref degree, ref minute, ref second);
betaOrigin = betaOrigin * DegToRad;
lambdaOrigin = lambdaOrigin * DegToRad;
try
{
// X,Yの入力
if (ellipSuffix == 1)
{
X = inLat;
Y = inLon;
}
else
{
X = inLat;
Y = inLon;
}
// 赤道からの子午線長の計算
MeridS(betaOrigin, AEE, AJ, BJ, CJ, DJ, EJ, FJ, GJ, HJ, IJ, ref S0); // 赤道から座標原点までの子午線弧長
M = S0 + X / M0;
// Baileyの式による異性緯度（isometric latitude）phi1の計算。
// 「精密測地網一次基準点測量計算式」P57の11.(1)の式から。
// この式と「現代測量学１ 測量の数学的基礎」P102の式とは，Math.Cos(phi1)だけ異なる。
// この式を導入したためベッセル楕円体以外で往復計算OKとなった。
icount = 0;
phi1 = betaOrigin;
do
{
icount = icount + 1;
oldphi1 = phi1;
MeridS(phi1, AEE, AJ, BJ, CJ, DJ, EJ, FJ, GJ, HJ, IJ, ref S1);
Bunsi = 2.0 * (S1 - M) * Math.Pow((1.0 - EPs.Eccentricity * Math.Sin(phi1) * Math.Sin(phi1)), 1.5);
Bunbo = 3.0 * EPs.Eccentricity * (S1 - M) * Math.Sin(phi1) * Math.Cos(phi1) * Math.Sqrt(1.0 - EPs.Eccentricity * Math.Sin(phi1) * Math.Sin(phi1)) - 2.0 * EPs.Axix * (1.0 - EPs.Eccentricity);
phi1 = phi1 + Bunsi / Bunbo;
if (icount > 100)
MessageBox.Show("異性緯度を求めるのに収束しません。", Application.ProductName, MessageBoxButtons.OK, MessageBoxIcon.Error);
}
while (!((Math.Abs(phi1 - oldphi1) < 0.00000000000001) || (icount > 100))); // 赤道から点までの子午線弧長
// 本では1e-12で十分 iterationの回数は４回
// 何度も使う式を変数に代入
double YM0 = Y / M0;
double T = Math.Tan(phi1); // 「精密測地網一次基準点測量計算式」P51のt1に等しい
double T2 = T * T; T4 = T2 * T2; T6 = T4 * T2;
double Eta2 = Ep2 * Math.Cos(phi1) * Math.Cos(phi1); // =η1*η1
double M1 = CEE / Math.Sqrt(Math.Pow((1.0 + Eta2), 3.0));
double N1 = CEE / Math.Sqrt(1.0 + Eta2);
double N1CosPhi1 = N1 * Math.Cos(phi1);
// 緯度Bの計算 「精密測地網一次基準点測量計算式」P51のphiを求める式より
beta = ((1385.0 + 3633.0 * T2 + 4095 * T4 + 1575.0 * T6) / (40320.0 * Math.Pow(N1, 8.0))) * Math.Pow(YM0, 8.0);
beta = beta - ((61.0 + 90.0 * T2 + 45 * T4 + 107.0 * Eta2 - 162.0 * T2 * Eta2 - 45.0 * T4 * Eta2) / (720.0 * Math.Pow(N1, 6.0))) * Math.Pow(YM0, 6.0);
beta = beta + ((5.0 + 3.0 * T2 + 6.0 * Eta2 - 6.0 * T2 * Eta2 - 3.0 * Math.Pow(Eta2, 2) - 9.0 * T2 * Math.Pow(Eta2, 2)) / (24.0 * Math.Pow(N1, 4.0))) * Math.Pow(YM0, 4.0);
beta = beta - ((1.0 + Eta2) / (2.0 * Math.Pow(N1, 2.0))) * Math.Pow(YM0, 2.0);
beta = beta * T;
beta = beta + phi1;
// 経度Lの計算 「精密測地網一次基準点測量計算式」P51のΔλを求める式より
lambda = -((61.0 + 662.0 * T2 + 1320.0 * T4 + 720.0 * T6) / (5040.0 * Math.Pow(N1, 6.0) * N1CosPhi1)) * Math.Pow(YM0, 7.0);
lambda = lambda + ((5.0 + 28.0 * T2 + 24.0 * T4 + 6.0 * Eta2 + 8.0 * T2 * Eta2) / (120.0 * Math.Pow(N1, 4.0) * N1CosPhi1)) * Math.Pow(YM0, 5.0);
lambda = lambda - ((1.0 + 2.0 * T2 + Eta2) / (6.0 * Math.Pow(N1, 2.0) * N1CosPhi1)) * Math.Pow(YM0, 3.0);
lambda = lambda + (1.0 / N1CosPhi1) * YM0;
lambda = lambda + lambdaOrigin;
// 子午線収差角の計算 「精密測地網一次基準点測量計算式」P51のγを求める式より
gamma = ((1.0 + T2) * (2.0 + 3.0 * T2) / (15.0 * Math.Pow(N1, 5.0))) * Math.Pow((Y / M0), 5.0);
gamma = gamma - ((1.0 + T2 - Eta2) / (3.0 * Math.Pow(N1, 3.0))) * Math.Pow((Y / M0), 3.0);
gamma = gamma + (1.0 / N1) * Y / M0;
gamma = gamma * T;
// 縮尺係数の計算 「精密測地網一次基準点測量計算式」P51のmを求める式より
double Eta2phi = Ep2 * Math.Cos(beta) * Math.Cos(beta); // =η*η。Bはphiと同じ。
double Mphi = CEE / Math.Sqrt(Math.Pow((1.0 + Eta2phi), 3.0));
double Nphi = CEE / Math.Sqrt(1.0 + Eta2phi);
double MMM = 0.0; // 縮尺係数
MMM = Math.Pow(Y, 4.0) / (24.0 * Mphi * Mphi * Nphi * Nphi * Math.Pow(M0, 4.0));
MMM = MMM + Y * Y / (2.0 * Mphi * Nphi * Math.Pow(M0, 2.0));
MMM = MMM + 1.0;
MMM = MMM * M0;
// 出力
double betaDeg = 0.0; // 求める緯度，経度。基本的にdeg
betaDeg = beta * RadToDeg;
double lambdaDeg = 0.0; // 求める緯度，経度。基本的にdeg
lambdaDeg = lambda * RadToDeg;
// 子午線収差角Gamma(rad)から，真北方向角Gammadeg(degree)へ
double gammaDeg = 0.0; // 真北方向角。deg
gammaDeg = -gamma * RadToDeg;
//緯度、経度
DegToDms(betaDeg, ref degree, ref minute, ref second, ref signString);
if (ellipSuffix == 1)
outLat = $"{signString}{degree:###}{minute:00}{second:00.00000}";
else
outLon = $"{signString}{degree:###}{minute:00}{second:00.00000}";
//X,Y
DegToDms(lambdaDeg, ref degree, ref minute, ref second, ref signString);
if (ellipSuffix == 1)
outX = $"{signString}{degree:###}{minute:00}{second:00.00000}";
else
outY = $"{signString}{degree:###}{minute:00}{second:00.00000}";
//真北方向角
DegToDmsNotZero3(gammaDeg, ref degree, ref minute, ref second, ref signString);
outDeg = $"{signString}{degree:###}{minute:00}{second:00.000}";
//縮尺係数
outParam = $"{MMM:#.00000000}"; // 成果表は小数点以下６桁
// outLat, outLon, outX, outY, outDeg, outParam に出力値が入っているので呼び出し下で受け取れるように調整を。
}
catch (OverflowException)
{
CommonUtils.LogWrite(string.Format(ConstString.Error.FailedXY, MethodBase.GetCurrentMethod().Name, beta, lambda));
return;
}
catch (Exception e)
{
CommonUtils.LogWrite(string.Format(ConstString.Error.FailedXY, MethodBase.GetCurrentMethod().Name, beta, lambda));
return;
}
}
/// <summary>
/// 赤道から緯度Phiまでの子午線弧長を計算します。
/// 「精密測地網一次基準点測量計算式」P55より
/// </summary>
/// <param name="Phi"></param>
/// <param name="AEE"></param>
/// <param name="AJ"></param>
/// <param name="BJ"></param>
/// <param name="CJ"></param>
/// <param name="DJ"></param>
/// <param name="EJ"></param>
/// <param name="FJ"></param>
/// <param name="GJ"></param>
/// <param name="HJ"></param>
/// <param name="IJ"></param>
/// <param name="SS"></param>
public void MeridS(double Phi, double AEE, double AJ, double BJ, double CJ, double DJ, double EJ, double FJ, double GJ, double HJ, double IJ, ref double SS)
{
SS = IJ / 16.0 * Math.Sin(16.0 * Phi);
SS = SS - HJ / 14.0 * Math.Sin(14.0 * Phi);
SS = SS + GJ / 12.0 * Math.Sin(12.0 * Phi);
SS = SS - FJ / 10.0 * Math.Sin(10.0 * Phi);
SS = SS + EJ / 8.0 * Math.Sin(8.0 * Phi);
SS = SS - DJ / 6.0 * Math.Sin(6.0 * Phi);
SS = SS + CJ / 4.0 * Math.Sin(4.0 * Phi);
SS = SS - BJ / 2.0 * Math.Sin(2.0 * Phi);
SS = SS + AJ * Phi;
SS = AEE * SS;
}
/// <summary>
/// 「3ﾊﾟﾗﾒｰﾀによる」Tokyo97系からITRF94系への座標変換を行う。
/// この変換では楕円体高Hはゼロとする。
/// </summary>
/// <param name="B1">入力 緯度(度)</param>
/// <param name="L1">入力 経度(度)</param>
/// <param name="B2">出力 緯度(度)</param>
/// <param name="L2">出力 経度(度)</param>
public void Tokyo97toITRF94(double B1, double L1, ref double B2, ref double L2)
{
double Brad = 0.0;
double ALrad = 0.0; // 緯度，経度(radian)
double H = 0.0; // 楕円体高(m)
double X1 = 0.0;
double Y1 = 0.0;
double Z1 = 0.0;
double X2 = 0.0;
double Y2 = 0.0;
double Z2 = 0.0; // (m)
Brad = B1 * DegToRad;
ALrad = L1 * DegToRad;
// 緯度，経度から三次元直交座標系(X,Y,Z)への換算
BLHXYZcalc(Brad, ALrad, 0, ref X1, ref Y1, ref Z1, ellipObj[0]); // EP(0):Bessel楕円体=日本測地系
// 三次元直交座標系(X,Y,Z)から三次元直交座標系(X,Y,Z)への座標変換
Xyz2Xyz(X1, Y1, Z1, ref X2, ref Y2, ref Z2);
// 三次元直交座標系(X,Y,Z)から緯度，経度への換算
XYZBLHcalc(X2, Y2, Z2, ref Brad, ref ALrad, ref H, ellipObj[1]); // EP(1):GRS-80楕円体
B2 = Brad * RadToDeg;
L2 = ALrad * RadToDeg;
}
/// <summary>
/// 「3ﾊﾟﾗﾒｰﾀによる」ITRF94系からTokyo97系への座標変換を行う。
/// この変換では楕円体高Hは，変換後の水平位置B2,L2を正確に元に戻すため，2ステップで変換する。
/// (1)初期値はゼロとし一度変換したあと，(2)その時のH分だけ下駄を履かせ，もう一度変換する。
/// </summary>
/// <param name="B1">入力 緯度(度)</param>
/// <param name="L1">入力 経度(度)</param>
/// <param name="B2">出力 緯度(度)</param>
/// <param name="L2">出力 経度(度)</param>
public void ITRF94toTokyo97(double B1, double L1, ref double B2, ref double L2)
{
double Brad1 = 0.0; // 緯度，経度(radian)
double ALrad1 = 0.0; // 緯度，経度(radian)
double Brad2 = 0.0; // 緯度，経度(radian)
double ALrad2 = 0.0; // 緯度，経度(radian)
double H = 0.0; // 楕円体高(m)
double X1 = 0.0;
double Y1 = 0.0;
double Z1 = 0.0;
double X2 = 0.0;
double Y2 = 0.0;
double Z2 = 0.0; // (m)
Brad1 = B1 * DegToRad;
ALrad1 = L1 * DegToRad;
// (1)
// 緯度，経度から三次元直交座標系(X,Y,Z)への換算
BLHXYZcalc(Brad1, ALrad1, 0.0, ref X1, ref Y1, ref Z1, ellipObj[1]); // EP(1):GRS-80楕円体=世界測地系
// 三次元直交座標系(X,Y,Z)から三次元直交座標系(X,Y,Z)への座標変換
Xyz2XyzR(X1, Y1, Z1, ref X2, ref Y2, ref Z2);
// 三次元直交座標系(X,Y,Z)から緯度，経度への換算
XYZBLHcalc(X2, Y2, Z2, ref Brad2, ref ALrad2, ref H, ellipObj[0]); // EP(0):Bessel楕円体=日本測地系
// (2)
// 緯度，経度から三次元直交座標系(X,Y,Z)への換算
BLHXYZcalc(Brad1, ALrad1, -H, ref X1, ref Y1, ref Z1, ellipObj[1]); // EP(0):GRS-80楕円体=世界測地系
// 三次元直交座標系(X,Y,Z)から三次元直交座標系(X,Y,Z)への座標変換
Xyz2XyzR(X1, Y1, Z1, ref X2, ref Y2, ref Z2);
// 三次元直交座標系(X,Y,Z)から緯度，経度への換算
XYZBLHcalc(X2, Y2, Z2, ref Brad2, ref ALrad2, ref H, ellipObj[0]); // EP(1):Bessel楕円体=日本測地系
B2 = Brad2 * RadToDeg;
L2 = ALrad2 * RadToDeg;
}
/// <summary>
/// Ver.1.2から，EllipPar Typeに対応
/// B,L,HからX,Y,Zに変換する。 このとき楕円体を指定する必要あり。
/// </summary>
/// <param name="Brad">緯度(RADIAN)</param>
/// <param name="ALrad">経度(RADIAN)</param>
/// <param name="H">高さ(m)</param>
/// <param name="X">３次元直交座標系(m)</param>
/// <param name="Y">３次元直交座標系(m)</param>
/// <param name="Z">３次元直交座標系(m)</param>
/// <param name="EP">計算に利用する楕円体定義</param>
private void BLHXYZcalc(double Brad, double ALrad, double H, ref double X, ref double Y, ref double Z, EllipParam EP) // ________________________
{
double CB = 0.0;
double CL = 0.0;
double SB = 0.0;
double SL = 0.0;
double W = 0.0;
double AN = 0.0;
double a = 0.0;
double FI = 0.0;
double e = 0.0;
a = EP.Axix;
FI = EP.Flattening1;
e = EP.Eccentricity;
CB = Math.Cos(Brad); SB = Math.Sin(Brad);
CL = Math.Cos(ALrad); SL = Math.Sin(ALrad);
W = Math.Sqrt(1.0 - e * SB * SB);
AN = a / W;
X = (AN + H) * CB * CL;
Y = (AN + H) * CB * SL;
Z = (AN * (1.0 - e) + H) * SB;
}
/// <summary>
/// Ver.2.2から，EllipPar Typeに対応
/// X,Y.ZからB,L,Hに変換する。 このとき楕円体を指定する必要あり。
/// </summary>
/// <param name="X">３次元直交座標系(m)</param>
/// <param name="Y">３次元直交座標系(m)</param>
/// <param name="Z">３次元直交座標系(m)</param>
/// <param name="Brad">緯度(RADIAN)</param>
/// <param name="ALrad">経度(RADIAN)</param>
/// <param name="H">高さ(m)</param>
/// <param name="EP">楕円体定数</param>
private void XYZBLHcalc(double X, double Y, double Z, ref double Brad, ref double ALrad, ref double H, EllipParam EP)
{
double p2 = 0.0;
double p = 0.0;
double r = 0.0;
double myu = 0.0;
double myus3 = 0.0;
double myuc3 = 0.0;
double a = 0.0;
double FI = 0.0;
double e = 0.0;
a = EP.Axix;
FI = EP.Flattening1;
e = EP.Eccentricity;
p2 = X * X + Y * Y;
p = Math.Sqrt(p2);
if (p == 0)
{
// エラー処理ルーチン。
MessageBox.Show(
"Warning: 座標値(X,Y,Z)を入力して下さい。",
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Error
);
return;
}
if (X == 0)
ALrad = HalfPI;
else
ALrad = Math.Atan(Y / X);// 経度
if (X < 0)
ALrad = ALrad + PI;
if (ALrad > PI)
ALrad = ALrad - DoublePI;
r = Math.Sqrt(p2 + Z * Z);
myu = Math.Atan((Z / p) * (1.0 - (1.0 / FI) + e * a / r));
myus3 = Math.Sin(myu);
myus3 = myus3 * myus3 * myus3;
myuc3 = Math.Cos(myu);
myuc3 = myuc3 * myuc3 * myuc3;
Brad = Math.Atan((Z * (1.0 - 1.0 / FI) + e * a * myus3) / ((1.0 - 1.0 / FI) * (p - e * a * myuc3)));
H = p * Math.Cos(Brad) + Z * Math.Sin(Brad) - a * Math.Sqrt(1.0 - e * Math.Sin(Brad) * Math.Sin(Brad)); // 楕円体高
}
/// <summary>
/// 座標変換をします。
/// 変換ﾊﾟﾗﾒｰﾀは，Tokyo97からITRF94固定とします。TKY2JGD専用なので。
/// |XS| |X| |T1| || D -R3 R2||X|
/// |YS|=|Y|+|T2|+| R3 D -R1||Y| cf. IERS ANNUAL REPORT FOR 1989
/// |ZS| |Z| |T3| |-R2 R1 D ||Z| II-23
/// </summary>
/// <param name="X1">３次元直交座標系での座標 (m)</param>
/// <param name="Y1">３次元直交座標系での座標 (m)</param>
/// <param name="Z1">３次元直交座標系での座標 (m)</param>
/// <param name="X2">３次元直交座標系での座標 (m)</param>
/// <param name="Y2">３次元直交座標系での座標 (m)</param>
/// <param name="Z2">３次元直交座標系での座標 (m)</param>
public void Xyz2Xyz(double X1, double Y1, double Z1, ref double X2, ref double Y2, ref double Z2)
{
const double T1 = -14641.4; // X shift(cm)
const double T2 = 50733.7; // Y shift(cm)
const double T3 = 68050.7; // Z shift(cm)
double T1real = 0.0; // X shift(m)
double T2real = 0.0; // Y shift(m)
double T3real = 0.0; // Z shift(m)
// パラメータの代入 Tokyo97からITRF94へ
// -14641.40 50733.70 68050.70 0.00 0.00 0.00 0.00 Tokyo97 ITRF94
// 実際に使う変換パラメータの計算
T1real = T1 * 0.01; // cm to m
T2real = T2 * 0.01; // cm to m
T3real = T3 * 0.01; // cm to m
// 一般的には７パラメータで次のように変換しますが，ここでは3パラメータなので不要な計算は省略します。
X2 = X1 + T1real;
Y2 = Y1 + T2real;
Z2 = Z1 + T3real;
}
/// <summary>
/// 座標変換をします。
/// 逆方向変換 R:Reverse
/// 変換ﾊﾟﾗﾒｰﾀは，ITRF94からTokyo97固定とします。TKY2JGD専用なので。
/// 座標変換をします。
/// |XS| |X| |T1| || D -R3 R2||X|
/// |YS|=|Y|+|T2|+| R3 D -R1||Y| cf. IERS ANNUAL REPORT FOR 1989
/// |ZS| |Z| |T3| |-R2 R1 D ||Z| II-23
/// </summary>
/// <param name="X1">３次元直交座標系での座標 (m)</param>
/// <param name="Y1">３次元直交座標系での座標 (m)</param>
/// <param name="Z1">３次元直交座標系での座標 (m)</param>
/// <param name="X2">３次元直交座標系での座標 (m)</param>
/// <param name="Y2">３次元直交座標系での座標 (m)</param>
/// <param name="Z2">３次元直交座標系での座標 (m)</param>
public void Xyz2XyzR(double X1, double Y1, double Z1, ref double X2, ref double Y2, ref double Z2)
{
// パラメータの代入 ITRF94からTokyo97へ
// -14641.40 50733.70 68050.70 0.00 0.00 0.00 0.00 Tokyo97 ITRF94
const double T1 = 14641.4; // X shift(cm)
const double T2 = -50733.7; // Y shift(cm)
const double T3 = -68050.7; // Z shift(cm)
double T1real = 0.0; // X shift(cm)
double T2real = 0.0; // Y shift(cm)
double T3real = 0.0; // Z shift(cm)
// 実際に使う変換パラメータの計算
T1real = T1 * 0.01; // cm to m
T2real = T2 * 0.01; // cm to m
T3real = T3 * 0.01; // cm to m
// 一般的には７パラメータで次のように変換しますが，ここでは3パラメータなので不要な計算は省略します。
X2 = X1 + T1real;
Y2 = Y1 + T2real;
Z2 = Z1 + T3real;
}
/// <summary>
/// 度単位の緯度lat,経度lonを受け取り，その点が含まれるﾒｯｼｭコードを返す関数
/// もし境界線上の点が与えられると，その点は点の北や東のﾒｯｼｭに属すると解釈する
/// Public Type MeshCode
/// </summary>
/// <param name="lat">緯度</param>
/// <param name="lon">経度</param>
/// <returns></returns>
private MeshCode LatLon2MeshCode(double lat, double lon)
{
int lat1 = 0; // 1次メッシュコード
int lon1 = 0; // 1次メッシュコード
int LatLon1 = 0; // 1次メッシュコード
int lat2 = 0; // 2次メッシュコード
int lon2 = 0; // 2次メッシュコード
int latlon2 = 0; // 2次メッシュコード
int lat3 = 0; // 3次メッシュコード
int lon3 = 0; // 3次メッシュコード
int latlon3 = 0; // 3次メッシュコード
double ModLat = 0.0;
double ModLon = 0.0;
MeshCode rtn_mesh = new MeshCode();
// 1次メッシュコード各2桁
lat1 = (int)(lat * 1.5); // 36.0 -> 54, 36.666 -> 54
lon1 = (int)(lon) - 100; // 138.0 -> 38, 138.999 -> 38
// 2次メッシュコード各1桁
lat2 = (int)(8.0 * (1.5 * lat - lat1)); // 36.0 -> 0, 36.0833 -> 0
lon2 = (int)(8.0 * (lon - (lon1 + 100))); // 138.0 -> 0, 138.1249 -> 0
// 3次メッシュコード各1桁
lat3 = (int)(10.0 * (12.0 * lat - 8.0 * lat1 - lat2) + 0.00000000001);
lon3 = (int)(10.0 * (8.0 * (lon - (lon1 + 100)) - lon2) + 0.00000000001);
// 上2行が微少量を加えるのは，lon=138.45のとき3次が5とならずに6となるように
// Ver.2.1 to Ver.2.2 lat=36.0833333333333のときlat3=10になって，エラー画面が出るバグを解消。原因は微小量を加えたことの副作用。
// 本来は微小量加算をやめるべきだが，以前の計算値との継続性を重視し，正しい緯(経)度同士の繰り上がり処理で対応する。2002/02/21
if (lat3 == 10)
{
lat2 = lat2 + 1; lat3 = 0;
if (lat2 == 8)
{
lat1 = lat1 + 1; lat2 = 0;
}
}
if (lon3 == 10)
{
lon2 = lon2 + 1; lon3 = 0;
if (lon2 == 8)
{
lon1 = lon1 + 1; lon2 = 0;
}
}
// 3次メッシの左下（南西角）点から何度ずれているか？
ModLat = 120.0 * lat - 80.0 * lat1 - 10 * lat2 - lat3; // 余り。3次メッシの左下（南西角）点からどれくらいずれているか。0(西端)以上1(東端)未満
ModLon = 80.0 * (lon - (lon1 + 100)) - 10 * lon2 - lon3; // 余り。3次メッシの左下（南西角）点からどれくらいずれているか。0(南端)以上1(北端)未満
// 計算値のチェック
if (lat2 < 0 || lat2 > 7)
MessageBox.Show(
"緯度の2次メッシュコードが0-7の範囲にありません。 " + lat2.ToString(),
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Error
);
if (lon2 < 0 || lon2 > 7)
MessageBox.Show(
"経度の2次メッシュコードが0-7の範囲にありません。 " + lon2.ToString(),
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Error
);
if (lat3 < 0 || lat3 > 9)
MessageBox.Show(
"緯度の3次メッシュコードが0-9の範囲にありません。 " + lat3.ToString(),
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Error
);
if (lon3 < 0 || lon3 > 9)
MessageBox.Show(
"経度の3次メッシュコードが0-9の範囲にありません。 " + lon3.ToString(),
mainForm.ProductName,
MessageBoxButtons.OK,
MessageBoxIcon.Error
);
// 各メッシュコードの桁位置を調整し3次メッシュコードを合成する
LatLon1 = lat1 * 100 + lon1; // ex. 5438
latlon2 = lat2 * 10 + lon2; // ex. 23 00 to 77
latlon3 = lat3 * 10 + lon3; // ex. 45 00 to 99
rtn_mesh.FirstMeshCode = LatLon1; // 5438
rtn_mesh.SecondMeshCode = latlon2; // 23
rtn_mesh.ThridMeshCode = latlon3; // 45
rtn_mesh.FullMeshCode = Convert.ToInt64(LatLon1) * 10000 + latlon2 * 100 + latlon3; // 54382345
rtn_mesh.FormedMeshCode = $"{LatLon1:0000}-{latlon2:00}-{latlon3:00}"; // "5438-23-45"
rtn_mesh.ModLat = ModLat;
rtn_mesh.ModLon = ModLon;
return rtn_mesh;
}
/// <summary>
/// 基本ﾒｯｼｭｺｰﾄﾞoriginMeshCodeを受け取り，隣の東，北，北東のﾒｯｼｭコードを返す[TonariMeshCode2]
/// Public Type MeshCode '参考のため構造体を以下に示す
/// called by Bilinear1 バイリニア補間用
/// </summary>
/// <param name="originMeshCode">原点のメッシュコード</param>
/// <param name="eastMeshCode">東隣メッシュコード</param>
/// <param name="northMeshCode">北隣メッシュコード</param>
/// <param name="northEastMeshCode">北東隣メッシュコード</param>
private void AdjacentMeshCode(MeshCode originMeshCode, ref MeshCode eastMeshCode, ref MeshCode northMeshCode, ref MeshCode northEastMeshCode)
{
MeshCode MCbuf = new MeshCode();
int lat1 = 0; // 1次メッシュコード
int lon1 = 0; // 1次メッシュコード
int lat2 = 0; // 2次メッシュコード
int lon2 = 0; // 2次メッシュコード
int lat3 = 0; // 3次メッシュコード
int lon3 = 0; // 3次メッシュコード
int lat1out = 0; // 1次メッシュコード
int lon1out = 0; // 1次メッシュコード
int latlon1out = 0; // 1次メッシュコード
int lat2out = 0; // 2次メッシュコード
int lon2out = 0; // 2次メッシュコード
int latlon2out = 0; // 2次メッシュコード
int lat3out = 0; // 3次メッシュコード
int lon3out = 0; // 3次メッシュコード
int latlon3out = 0; // 3次メッシュコード
// ①まず基本ﾒｯｼｭｺｰﾄﾞを分解し要素を取り出す。
// 1次メッシュコード各2桁 'ex.54382345
lat1 = originMeshCode.FirstMeshCode / 100; // ex.54
lon1 = originMeshCode.FirstMeshCode % 100; // ex.38
// 2次メッシュコード各1桁
lat2 = originMeshCode.SecondMeshCode / 10; // ex.2
lon2 = originMeshCode.SecondMeshCode % 10; // ex.3
// 3次メッシュコード各1桁
lat3 = originMeshCode.ThridMeshCode / 10; // ex.4
lon3 = originMeshCode.ThridMeshCode % 10; // ex.5
// 東側ﾒｯｼｭ
// ②出力変数に代入したあと，となりのﾒｯｼｭｺｰﾄﾞを計算する。
lat1out = lat1;
lon1out = lon1;
lat2out = lat2;
lon2out = lon2;
lat3out = lat3;
lon3out = lon3;
if (lon3 != 9)
lon3out = lon3 + 1;
else
{
lon3out = 0;
if (lon2 != 7)
lon2out = lon2 + 1;
else
{
lon2out = 0;
lon1out = lon1 + 1;
}
}
// ③各メッシュコードの桁位置を調整し3次メッシュコードを合成する
latlon1out = lat1out * 100 + lon1out; // ex. 5438
latlon2out = lat2out * 10 + lon2out; // ex. 23 00 to 77
latlon3out = lat3out * 10 + lon3out; // ex. 45 00 to 99
// ④結果を仮の変数に代入
MCbuf.FirstMeshCode = latlon1out; // 5438
MCbuf.SecondMeshCode = latlon2out; // 23
MCbuf.ThridMeshCode = latlon3out; // 45
MCbuf.FullMeshCode = Convert.ToInt64(latlon1out) * 10000 + latlon2out * 100 + latlon3out; // 54382345
MCbuf.FormedMeshCode = $"{latlon1out:0000}-{latlon2out:00}-{latlon3out:00}"; // "5438-23-45"
// ⑤仮の変数から本当の変数に代入
eastMeshCode = MCbuf;
// 北側ﾒｯｼｭ
// ②出力変数に代入したあと，となりのﾒｯｼｭｺｰﾄﾞを計算する。
lat1out = lat1;
lon1out = lon1;
lat2out = lat2;
lon2out = lon2;
lat3out = lat3;
lon3out = lon3;
if (lat3 != 9)
lat3out = lat3 + 1;
else
{
lat3out = 0;
if (lat2 != 7)
lat2out = lat2 + 1;
else
{
lat2out = 0;
lat1out = lat1 + 1;
}
}
// ③各メッシュコードの桁位置を調整し3次メッシュコードを合成する
latlon1out = lat1out * 100 + lon1out; // ex. 5438
latlon2out = lat2out * 10 + lon2out; // ex. 23 00 to 77
latlon3out = lat3out * 10 + lon3out; // ex. 45 00 to 99
// ④結果を仮の変数に代入
MCbuf.FirstMeshCode = latlon1out; // 5438
MCbuf.SecondMeshCode = latlon2out; // 23
MCbuf.ThridMeshCode = latlon3out; // 45
MCbuf.FullMeshCode = Convert.ToInt64(latlon1out) * 10000 + latlon2out * 100 + latlon3out; // 54382345
MCbuf.FormedMeshCode = $"{latlon1out:0000}-{latlon2out:00}-{latlon3out:00}"; // "5438-23-45"
// ⑤仮の変数から本当の変数に代入
northMeshCode = MCbuf;
// 北東側ﾒｯｼｭ
// ②出力変数に代入したあと，となりのﾒｯｼｭｺｰﾄﾞを計算する。
lat1out = lat1;
lon1out = lon1;
lat2out = lat2;
lon2out = lon2;
lat3out = lat3;
lon3out = lon3;
if (lon3 != 9)
lon3out = lon3 + 1;
else
{
lon3out = 0;
if (lon2 != 7)
lon2out = lon2 + 1;
else
{
lon2out = 0;
lon1out = lon1 + 1;
}
}
if (lat3 != 9)
lat3out = lat3 + 1;
else
{
lat3out = 0;
if (lat2 != 7)
lat2out = lat2 + 1;
else
{
lat2out = 0;
lat1out = lat1 + 1;
}
}
// ③各メッシュコードの桁位置を調整し3次メッシュコードを合成する
latlon1out = lat1out * 100 + lon1out; // ex. 5438
latlon2out = lat2out * 10 + lon2out; // ex. 23 00 to 77
latlon3out = lat3out * 10 + lon3out; // ex. 45 00 to 99
// ④結果を仮の変数に代入
MCbuf.FirstMeshCode = latlon1out; // 5438
MCbuf.SecondMeshCode = latlon2out; // 23
MCbuf.ThridMeshCode = latlon3out; // 45
MCbuf.FullMeshCode = Convert.ToInt64(latlon1out) * 10000 + latlon2out * 100 + latlon3out; // 54382345
MCbuf.FormedMeshCode = string.Format(latlon1out.ToString(), "0000") + "-" + string.Format(latlon2out.ToString(), "00") + "-" + string.Format(latlon3out.ToString(), "00"); // "5438-23-45"
// ⑤仮の変数から本当の変数に代入
northEastMeshCode = MCbuf;
}
}
}
使い方
Tky2jgd.Tky2jgd zahyou_converter = new Tky2jgd.Tky2jgd(stParentName01, this);
zahyou_converter.FileConvertBatch();
このようにやることで、 stParentName01 フォルダ内の、 zahyou.in ファイルを読み込み、 zahyou.out ファイルを変換結果として出力します。
第3～5引数を指定することで、zahyou.in zahyou.out TKY2JGD.par の各ファイル名を指定することもできます。
最後に
大元のソフトウェアを開発してくださった国土地理院の 飛田幹男 氏に感謝を。
またソフトウェアの改変については、国土地理院の規約 を確認の上、問題ないと判断し私が対応したものとなります。
C#版につきまして疑問が生じた場合には 国土地理院に対して質問等を行うことが無いようお願いいたします 。