- 投稿日：2020-03-24T20:14:15+09:00
警告なしでenumをforeach文で回す
- 投稿日：2020-03-24T18:00:15+09:00
C#でビンゴゲームのシミュレーション
C#でビンゴゲームのシミュレーションを行いました
ランダムに数字を選びビンゴカードを5万枚作り、ビンゴゲームをシミュレーションしました。
ビンゴカードは左の列から、B列、I列、N列、G列、O列といい、それぞれ数字の範囲が
1～15、16～30、31～45、46～60、61～75と決まっています。
カードの数字の配列を2次元配列で表しました。ビンゴマシーンクラス
名称は分かりませんが、ハンドル回してガラガラガラして、数字の書かれたボールを出す、あのマシーン。
それを一つのクラスにしましたので、ライブラリー化することもできますし、コピーして、プロジェクトの一部にもできます。BingoMachine.csusing System.Collections.Generic; namespace CsBingoTest { class BingoMachine { private List<int> number_list = new List<int>(75); public BingoMachine(System.Random rnd) //コンストラクターでRandomを引数にするのは、短時間に複数のRandomインスタンスを生成すると、Next()実行時に同じ数字が返されるため。 //従って、呼び出し元であらかじめ用意したRandomインスタンスを引数に渡して、複数のBingoMachineインスタンスを短時間に生成する場合は、共通のRandomインスタンスを引数に渡してください。 { var list = new List<int>(75); int l; for(int i = 1; i <= 75; i++) { list.Add(i); } for (int i = 0; i < 75; i++) { l = rnd.Next(0, 75 - i); number_list.Add(list[l]); list.RemoveRange(l, 1); } } public int rolling() //ビンゴを、ハンドル回してガラガラガラ～。出た数字を返します。 { if (number_list.Count <= 0) { return 0; } int num = number_list[0]; number_list.RemoveRange(0, 1); return num; } } }ビンゴカードクラス
コードは以下の通りです。
BINGO.csnamespace CsBingoTest { class BINGO { /// <summary> /// ビンゴカードの数字 /// </summary> private int[,] num = new int[5, 5] { { 0, 0, 0, 0, 0 }, { 0, 0, 0, 0, 0 }, { 0, 0, 0, 0, 0 }, { 0, 0, 0, 0, 0 }, { 0, 0, 0, 0, 0 } }; /// <summary> /// ビンゴカードに穴が開いたかどうかを記録 /// </summary> private bool[,] b = new bool[5, 5] { { false, false, false, false, false }, { false, false, false, false, false }, { false, false, true, false, false }, { false, false, false, false, false }, { false, false, false, false, false } }; public BINGO(System.Random rnd) //コンストラクターでRandomを引数にするのは、短時間に複数のRandomインスタンスを生成すると、Next()実行時に同じ数字が返されるため。 //従って、呼び出し元であらかじめ用意したRandomインスタンスを引数に渡して、複数のBINGOインスタンスを短時間に生成する場合は、共通のRandomインスタンスを引数に渡してください。 { int i; bool c; int num_p; int num_bingo_p; for (num_bingo_p = 0; num_bingo_p < 5; num_bingo_p++) { num_p = 0; while (num_p < 5) { c = false; i = rnd.Next(1, 15); for (int j = 0; j < 5; j++) { if (num[num_bingo_p, j] == i + (15 * num_bingo_p)) { c = true; } } if (!c) { num[num_bingo_p, num_p] = i + (15 * num_bingo_p); num_p++; } } } } public void push(int number) //ビンゴカードに穴を開ける関数です。intを渡すと、その数値があれば、穴を開けます。 { int bingo = 0; if (1 <= number && number <= 15) bingo = 0; else if (16 <= number && number <= 30) bingo = 1; else if (31 <= number && number <= 45) bingo = 2; else if (46 <= number && number <= 60) bingo = 3; else if (61 <= number && number <= 75) bingo = 4; for(int i = 0; i < 5; i++) { if (num[bingo, i] == number) { b[bingo, i] = true; } } } public bool isBingo //ビンゴしたかどうかを返します。 { get { bool bng = false; for(int i = 0; i < 5; i++) { if (b[i, 0] && b[i, 1] && b[i, 2] && b[i, 3] && b[i, 4]) { bng = true; } if (b[0, i] && b[1, i] && b[2, i] && b[3, i] && b[4, i]) { bng = true; } } if (b[0, 0] && b[1, 1] && b[2, 2] && b[3, 3] && b[4, 4]) { bng = true; } if (b[0, 4] && b[1, 3] && b[2, 2] && b[3, 1] && b[4, 0]) { bng = true; } return bng; } } public bool dummy = false; //何かのフラグにしてください。今回は、ビンゴした後にisBingoプロパティを参照したかどうかのフラグに使用してます。 } }メインプログラム
Program.csnamespace CsBingoTest { class Program { static void Main(string[] args) { int bingo_card = 50000; //ビンゴカードの枚数 System.Random rnd = new System.Random(); BingoMachine bingoMachine = new BingoMachine(rnd); BINGO[] bingo = new BINGO[bingo_card]; for(int i = 0; i < bingo_card; i++) { bingo[i] = new BINGO(rnd); //ビンゴカードの枚数分newする } int n; int loop = 1; //繰返し回数 int bingoed_num = 0; //1回あたりのビンゴした数 bool isAllBingo = false; while (!isAllBingo) //全てのカードがビンゴするまで繰り返す { isAllBingo = true; n = bingoMachine.rolling(); System.Console.Write(loop.ToString("00") + "=>" + n.ToString("00") + "::"); for (int i = 0; i < bingo_card; i++) { if (!bingo[i].isBingo) bingo[i].push(n); if (!bingo[i].isBingo) isAllBingo = false; if (bingo[i].isBingo && (!bingo[i].dummy)) { System.Console.Write(i.ToString() + ","); bingo[i].dummy = true; bingoed_num++; } } System.Console.Write("\n(" + bingoed_num.ToString() + ")\n"); bingoed_num = 0; loop++; } System.Console.WriteLine("ALL BINGO!!!"); } } }実行結果
実行結果は
回数=>出た数字::ビンゴしたカード番号...
(ビンゴしたカードの数)
となります。
ビンゴしたカードの番号も含めて全部ここに載せられないので、それを省いた実行結果を載せます。01=>39:: (0) 02=>44:: (0) 03=>59:: (0) 04=>28:: (0) 05=>02:: (0) 06=>19:: (0) 07=>05:: (0) 08=>56:: (0) 09=>70:: (38) 10=>04:: (13) 11=>63:: (59) 12=>07:: (30) 13=>60:: (0) 14=>42:: (15) 15=>48:: (79) 16=>20:: (106) 17=>33:: (81) 18=>32:: (210) 19=>13:: (123) 20=>61:: (270) 21=>08:: (271) 22=>36:: (526) 23=>38:: (1056) 24=>58:: (334) 25=>51:: (315) 26=>27:: (544) 27=>35:: (1796) 28=>62:: (747) 29=>46:: (664) 30=>37:: (2628) 31=>50:: (748) 32=>53:: (1150) 33=>49:: (1806) 34=>21:: (1116) 35=>30:: (0) 36=>55:: (2847) 37=>31:: (3270) 38=>41:: (4480) 39=>73:: (1019) 40=>11:: (980) 41=>15:: (0) 42=>03:: (1124) 43=>66:: (1249) 44=>09:: (1419) 45=>26:: (1457) 46=>18:: (1408) 47=>74:: (1454) 48=>14:: (1795) 49=>65:: (1381) 50=>10:: (1928) 51=>67:: (1173) 52=>43:: (2194) 53=>52:: (1135) 54=>24:: (843) 55=>71:: (723) 56=>06:: (998) 57=>12:: (1103) 58=>01:: (1295) ALL BINGO!!!5万枚あれば4，5回でビンゴするカードがあるかと思いましたが、この結果だけを見ると、9回目で初ビンゴしてます。
5万枚では、試行回数が少ないんでしょうか・・・？100万枚でやってみました。
01=>41:: (0) 02=>29:: (0) 03=>57:: (0) 04=>72:: (0) 05=>34:: (0) 06=>13:: (86) 07=>62:: (101) 08=>07:: (186) 09=>71:: (197) 10=>01:: (261) 11=>08:: (287) 12=>31:: (100) 13=>54:: (1293) 14=>49:: (1259) 15=>36:: (1314) 16=>47:: (1264) 17=>64:: (1706) 18=>70:: (2152) 19=>60:: (0) 20=>45:: (0) 21=>56:: (2578) 22=>26:: (10563) 23=>63:: (6695) 24=>46:: (7298) 25=>43:: (5735) 26=>48:: (12245) 27=>44:: (11602) 28=>59:: (21239) 29=>03:: (10843) 30=>38:: (20937) 31=>39:: (34767) 32=>18:: (26361) 33=>16:: (25754) 34=>52:: (39070) 35=>33:: (49750) 36=>30:: (0) 37=>42:: (70401) 38=>58:: (50462) 39=>06:: (20437) 40=>11:: (21605) 41=>32:: (85085) 42=>75:: (0) 43=>51:: (61492) 44=>02:: (20007) 45=>21:: (26974) 46=>04:: (25796) 47=>68:: (23087) 48=>15:: (0) 49=>67:: (23771) 50=>10:: (31118) 51=>22:: (29049) 52=>17:: (26591) 53=>25:: (24902) 54=>66:: (23739) 55=>28:: (22290) 56=>12:: (24574) 57=>40:: (25629) 58=>09:: (19526) 59=>24:: (10119) 60=>65:: (8662) 61=>37:: (13365) 62=>05:: (7450) 63=>14:: (8226) ALL BINGO!!!6回目で86枚ビンゴしてます。非常に低い確率だと分かります。
5回目でビンゴする確率はもっともっと低いのでしょう。
4回目でのビンゴは、もっともっともっともっともっと確率が低いのかもしれません。
- 投稿日：2020-03-24T16:28:27+09:00
【uGUI】 UI要素のサイズを別UI要素と一致させるスクリプト
背景
UI要素のサイズ(RectTransform.sizeDelta)を別UI要素に追従させ、見た目上で一致させたい
- 一致させる対象が直接の親の場合
- Stretch指定でOK
- 一致させる対象が子かつ子の数が1つの場合
- 下記公式ドキュメントの実装で可能
- https://docs.unity3d.com/ja/2018.4/Manual/HOWTO-UIFitContentSize.html
- Horizontal/Vertical Layout Groupを使用する必要があり、自動で整列されるため子1つのみの制約あり
- 構成が複雑でわかりづらい
上記の様なヒエラルキーや要素数に依存しないコンポーネントが欲しかったので自作しました
また、Inspector上でのサイズ変更を不可にしましたUISizeFitter
using UnityEngine; [ExecuteInEditMode, RequireComponent(typeof(RectTransform))] public class UISizeFitter : MonoBehaviour { [SerializeField] RectTransform _target; RectTransform _selfRect; Vector3 _cachedTargetLossyScale; Vector2 _cachedTargetSizeDelta; bool _isDirty = true; void UpdateDirty() { if (_cachedTargetSizeDelta == _target.sizeDelta && _cachedTargetLossyScale == _target.lossyScale) { return; } _cachedTargetSizeDelta = _target.sizeDelta; _cachedTargetLossyScale = _target.lossyScale; _isDirty = true; } void AdjustSize() { if (!_isDirty) { return; } _selfRect.sizeDelta = _cachedTargetSizeDelta * _cachedTargetLossyScale; _isDirty = false; } // Layout Groupみたいに「Some values driven by 〜」とInspector上の編集を禁止する void LockRectTransformProperty() { DrivenRectTransformTracker dt = new DrivenRectTransformTracker(); dt.Clear(); dt.Add(this, _selfRect, DrivenTransformProperties.SizeDelta); } void Awake() { _selfRect = GetComponent<RectTransform>(); LockRectTransformProperty(); } void Update() { if (_target == null) { return; } UpdateDirty(); AdjustSize(); } }補足
常に追従させたかったのでUpdateで実行していますが、負荷的に問題になったり、動的に変更しないのであればStart辺りに移動し一度のみ実行するのが良いと思います
また、lossyScaleがUnity側でキャッシュしているのか毎回計算しているのかが不明だったため、とりあえずComponent側でキャッシュしていますがUnity側で毎回計算している様だと比較時の負荷が若干気になります
_cached〜でキャッシュしているのは比較の方がコピーより高速だろうという考えの下です
- 投稿日：2020-03-24T15:12:04+09:00
今日からできる!わかりやすいプログラムを書く方法(C#)
概要
個人的にC#のプログラミングで気を付けていることです。一部は言語を問わないと思います。知っていれば、今日から簡単に気をつけることが可能なことを記載しています。
Windows Formsアプリを作成しているので偏っているかもしれません。個人的な好みもあると思いますし、正しいと主張しているわけではないので変なところがある場合にご指摘していただけると助かります。
1. 副作用を少なくする
個人的に一番主張したいことはこれです。
これは想定外の動きを少なくするために心掛けていることです。
自分の組んだプログラムでも2、3ヵ月たてば具体的な中身は忘れてしまいます。後々の苦労を減らすためにも副作用や影響範囲が明示的なプログラムを組むように心がけています。
副作用 (プログラム)1.1 staticな関数にする
staticな関数にすることにより、メンバー変数が思いがけず書き換えられることがないことを保証できます。極端ではありますがアクセスレベルがprivateでメンバー変数を触る必要のない関数は全部staticにしてしまえと思うぐらいです。
sample.cspublic int AnyProcess() { // メンバー変数が変更されるかどうかわからない } public static int AnyProcess() { // メンバー変数が変更されることはない }1.2 IReadOnlyListを使用する
中身が変更されることのないリストです。
IReadOnlyList インターフェイス引数をIReadOnlyListとして受け取ることで関数内でリストが変更されることないことがわかります。
sample.cspublic List<int> AnyProcess(List<int> list) { // 引数のListが変更されるかどうかがわからない } public List<int> AnyProcess(IReadOnlyList<int> list) { // 引数のListが変更されないことが明示的 }またプロパティの型をIReadOnlyListにすることにより、プロパティを参照する側は、このリストに対して直接AddしたりRemoveすることが想定されていないことがわかります。
sample.csprivate List<int> m_IntList; // そのまま返す場合 // 中身の変更が可能なので変更されてもいい場合以外はだめ public List<int> IntList { get { return m_IntList} } // Arrayに変換する場合 // 基本的にはこれで問題ないと思う public int[] IntList { get { return m_IntList.ToArray(); } } // IReadOnlyListとして公開する // 変更されたくないのが明示的。クラスの外からは基本的に変更されないはず // キャストすれば変更可能だけど...そこまでして変更する必要のある場合ともとれる？ public IReadOnlyList<int> IntList { get { return m_IntList; } }2. 公開範囲は小さくする
これはプログラムの変更をしやすくするための心得です。
後々変更するときに影響範囲を調べる必要が減ります。また変数が変更される範囲がわかりやすくなるので、どこで変更されているか調べる手間が減ります。2.1 必要最低限のアクセスレベルにする
関数や変数はprivateでいいならprivateで宣言して、internal⇒publicと
必要に応じて広げるようにします。
protectedは特殊ですが、意外と公開範囲は広いので気を付けてください
C#のアクセス修飾子 2019 〜protectedは 結構でかい〜2.2 private set
プロパティを公開したいことはあると思いますが、外からsetする必要がなければプロパティのsetterはprivate setにします。
sample.cs// クラスの外から変更可能 public int AnyValue { set; get; } // クラスの外から変更不可 public int AnyValue { private set; get; }2.3 readonly宣言とget-only property
コンストラクタ以外で変更する必要のない値はreadonlyをつけることにより、インスタンスが変更されないことが保証されます。（※Listの中身が変更されたり、クラスのメンバが変更されることはある）
また値を公開したい場合等はgetterだけの自動実装プロパティを用いることでコンストラクタ以外で変更されないことが保証できます。sample.cspublic class GetterSample { public double Name{ get; } public GetterSample(string name) { // 変更の必要がない場合get-onlyプロパティ Name= name; } }自動実装のget-onlyプロパティとpublic readonlyな変数の大きな違いはinterfaceで使用できるかどうかかなと思います。
sample.cspublic interface ISample { // プロパティはインターフェースでももてる string Hoge1 { get; } // これはだめ。フィールドはもてない readonly string Hoge2; }3. 見やすくする
見やすいってなんだって主観が入ることかもしれませんが、技術的にはそれほど難しくないことだけ。
3.1 適切にコメントをいれる
プログラムにコメントを書こうって話です。コメントの必要ないプログラムが書ければ、その方が良いと思いますが、一朝一夕ではいかないので、上手になるまではコメント書いた方がいいと思います。
コメントを書く基準を自分で決めると忘れにくいです。どういうときに書いたらいいかわからんって場合は、ググって誰かの参考にするか、リーダブルコード読むか、とりあえずコメントを書きまくって、しばらくして読み返していらないなってとこを消していって自分ルールを作っていけばいいと思います。
個人的には下記の2つは特に心がけてコメントするように心がけています。
・普通に考えたらこうするよねって記述ではなく一風変わった記述をしたときになぜそうしたかの理由・苦労して作って計算等、あとから見ると何をしているのかわからないような処理の説明
3.2 長い処理を（なるべく）書かない
"長くても100行まで。それ以上は処理を分ける"というように自分が見やすいルールを決めて書くといいと思います。チームで明確な基準があればいいと思うのですが、担当してる箇所によってはどうしても長くなってしまうようなケースもあると思うので基準がない場合は自分の見やすいようにすればいいと思う。
まとめ
ほかにもDRYとか変数や関数の命名とかいろいろ気を付けることはあると思いますが、簡単にできるようになるのは難しいのでここには記載していません。（特に変数や関数の命名ほんと難しいですよね。）
仕様変更に強く、使いやすい処理を書いて幸福なプログラミングライフを送りたい所存です。
- 投稿日：2020-03-24T14:12:40+09:00
日本測地系⇔世界測地系座標変換ソフト TKY2JGD をC#に
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
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.inzahyou.outTKY2JGD.parの各ファイル名を指定することもできます。最後に
大元のソフトウェアを開発してくださった国土地理院の 飛田幹男 氏に感謝を。
またソフトウェアの改変については、国土地理院の規約 を確認の上、問題ないと判断し私が対応したものとなります。
C#版につきまして疑問が生じた場合には 国土地理院に対して質問等を行うことが無いようお願いいたします 。
- 投稿日：2020-03-24T11:33:42+09:00
JupyterにMicrosoft.Data.Analysisを入れてみる
Jupyter環境に展開したC#にMicorosoft.Data.Analysisを載せてみました。
#r "nuget:Microsoft.Data.Analysis, 0.2.0"Installing package Microsoft.Data.Analysis, version 0.2.0.........done!
using Microsoft.Data.Analysis;#r "nuget:MathNet.Numerics,4.9.0"Installing package MathNet.Numerics, version 4.9.0.......done!
Successfully added reference to package MathNet.Numerics, version#r "nuget:XPlot.Plotly"Installing package XPlot.Plotly..........done!
Successfully added reference to package XPlot.Plotly, versionusing MathNet.Numerics; using XPlot.Plotly;var df = Microsoft.Data.Analysis.DataFrame.LoadCsv("Data/german.csv"); dfColumns Rows
[ [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ... (990 more) ], [ A11, A12, A14, A11, A11, A14, A14, A12, A14, A12 ... (990 more) ], [ 6, 48, 12, 42, 24, 36, 24, 36, 12, 30 ... (990 more) ], [ A34, A32, A34, A32, A33, A32, A32, A32, A32, A34 ... (990 more) ], [ A43, A43, A46, A42, A40, A46, A42, A41, A43, A40 ... (990 more) ], [ 1169, 5951, 2096, 7882, 4870, 9055, 2835, 6948, 3059, 5234 ... (990 more) ], [ A65, A61, A61, A61, A61, A65, A63, A61, A64, A61 ... (990 more) ], [ A75, A73, A74, A74, A73, A73, A75, A73, A74, A71 ... (990 more) ], [ 4, 2, 2, 2, 3, 2, 3, 2, 2, 4 ... (990 more) ], [ A93, A92, A93, A93, A93, A93, A93, A93, A91, A94 ... (990 more) ] ... (12 more) ] [ [ 0, A11, 6, A34, A43, 1169, A65, A75, 4, A93 ... (12 more) ], [ 1, A12, 48, A32, A43, 5951, A61, A73, 2, A92 ... (12 more) ], [ 2, A14, 12, A34, A46, 2096, A61, A74, 2, A93 ... (12 more) ], [ 3, A11, 42, A32, A42, 7882, A61, A74, 2, A93 ... (12 more) ], [ 4, A11, 24, A33, A40, 4870, A61, A73, 3, A93 ... (12 more) ], [ 5, A14, 36, A32, A46, 9055, A65, A73, 2, A93 ... (12 more) ], [ 6, A14, 24, A32, A42, 2835, A63, A75, 3, A93 ... (12 more) ], [ 7, A12, 36, A32, A41, 6948, A61, A73, 2, A93 ... (12 more) ], [ 8, A14, 12, A32, A43, 3059, A64, A74, 2, A91 ... (12 more) ], [ 9, A12, 30, A34, A40, 5234, A61, A71, 4, A94 ... (12 more) ] ... (990 more) ]となります。
さらに、可視化してみます。
var chart = Chart.Plot(new Graph.Scatter() {y = df["Credit amount"]});chartサンプルのJupyter Notebookをgithubに置いておきます。
- 投稿日：2020-03-24T10:07:10+09:00
WPF にしかない便利機能、階層を持ったデータを表示する HierarchicalDataTemplate
久しぶりの WPF ネタです。WPF から始まった MS 系の XAML で画面を定義するテクノロジーの中で、後発のプラットフォームでは採用されていない機能がいくつかあります。
MultiBinding とか、今回紹介する HierarchicalDataTemplate がそれにあたります。今回は HierarchicalDataTemplate に注目していきます。
前提
WPF では、データを画面に表示するときに、どのように表示するかという定義を DataTemplate というものを使って書きます。例えば以下のような Person クラスがあったとします。
namespace Xxx { public class Person { public string Name { get; set; } } }これを画面に文字列として Name を表示したいと思ったら以下のような感じのテンプレートになります。
<!-- 何処かで xmlns:xxx="clr-namespace:Xxx" を定義している前提 --> <DataTemplate TargetType="xxx:Person"> <TextBlock Text="{Binding Name}" /> </DataTemplate>DataTemplate 対応コントロールには ContentTemplate や ItemTemplate といったプロパティがあって、ここに上記の DataTemplate の定義を設定することで見た目を指定できます。
HierarchicalDataTemplate
HierarchicalDataTemplate は、階層構造を持つようなデータの見た目を定義するためのテンプレートになります。例えば以下のような Person クラスがあったとします。
using System.Collections.ObjectModel; namespace WpfApp40 { public class Person { public string Name { get; set; } public ObservableCollection<Person> Children { get; } = new ObservableCollection<Person>(); } }先ほどの例と比較すると Children プロパティが増えていますね！！Name を表示したあとに子要素として Children の Name も列挙したいというようなケースを考えてみましょう。階層構造のデータの表示に対応したコントロールとしては TreeView や ContextMenu コントロールがあります。今回は諸事情により ContextMenu でやってみようと思います。
MainPage.xaml.cs の DataContext に先ほどのクラスの配列を入れておきます。
using System.Windows; namespace WpfApp40 { public partial class MainWindow : Window { public MainWindow() { InitializeComponent(); DataContext = new[] { new Person { Name = "Taro1", Children = { new Person { Name = "Taro1-1", Children = { new Person { Name = "Taro1-1-1" }, new Person { Name = "Taro1-1-2" }, new Person { Name = "Taro1-1-3" }, new Person { Name = "Taro1-1-4" }, } }, new Person { Name = "Taro1-2" }, new Person { Name = "Taro1-3" }, }, }, new Person { Name = "Taro2", Children = { new Person { Name = "Taro2-1" }, new Person { Name = "Taro2-2" }, new Person { Name = "Taro2-3" }, }, }, }; } } }そして、画面にボタンを置いて、そのボタンのコンテキストメニューに HierarchicalDataTemplate を指定します。ポイントは DataTemplate の時と同様に Name プロパティを TextBlock で表示しているのに加えて ItemsSource プロパティで子要素として Children プロパティを指定しているところです。こうすることで、子要素には Children の中の Person クラスが列挙されて、その見た目が再度 HierarchicalDataTemplate で表示されるといった形になります。
<Window x:Class="WpfApp40.MainWindow" xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation" xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml" xmlns:d="http://schemas.microsoft.com/expression/blend/2008" xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006" xmlns:local="clr-namespace:WpfApp40" mc:Ignorable="d" Title="MainWindow" Height="450" Width="800"> <Grid HorizontalAlignment="Center" VerticalAlignment="Center"> <Button Content="XXXX" MinWidth="250"> <Button.ContextMenu> <ContextMenu ItemsSource="{Binding}"> <ContextMenu.ItemTemplate> <HierarchicalDataTemplate ItemsSource="{Binding Children}"> <TextBlock Text="{Binding Name}" /> </HierarchicalDataTemplate> </ContextMenu.ItemTemplate> </ContextMenu> </Button.ContextMenu> </Button> </Grid> </Window>実行すると以下のようになります。
これに DataTemplateSelector とかを組み合わせると、割と自由自在に階層構造を持ったデータを表現できて捗ります。
まとめ
WPF 強い。
- 投稿日：2020-03-24T08:51:25+09:00
Jupyter NotebookにC#をインストールする
.NET Core with Jupyter Notebooksにある通り、.NET CoreをJupyter Notebookから使えるようになったので、手元のWindows 環境でやってみました。Jupyter Notebookのインストールそのものは、Ububtu 18.04 環境上にJupyterとQ#環境を構築するあたりとか、その元とかが使えると思います。
基本的には
dotnet tool install -g --add-source "https://dotnet.myget.org/F/dotnet-try/api/v3/index.json" Microsoft.dotnet-interactive
とやって、.NET Interactiveをインストールして、Windows Terminalを再起動して、サーチパスを有効化し、
dotnet interactive jupyter installでjupyterのKernelをインストールしただけです。
あとは、おもむろにjupyter notebookでJupyter Notebookを起動。
System.Console.WriteLine("Foo");
とやって、Fooが出ることを確認できました。
うむ、すごく簡単。
次は、DataFrameで遊んでみよう。
- 投稿日：2020-03-24T06:25:19+09:00
つくるオーオース Introspect編
Client Credentials Grantと、Resource Owner Password Credentials Grantが再現出来たので、Introspectも再現しようと思います。
https://qiita.com/namikitakeo/items/38be899790cb27a323dfIntrospectコントローラーを作成します。
http://localhost:5000/op/introspectControllers/IntrospectController.csusing System; using System.IO; using System.Collections.Generic; using System.Linq; using System.Threading.Tasks; using Microsoft.AspNetCore.Http; using Microsoft.AspNetCore.Mvc; using Microsoft.EntityFrameworkCore; using myop.Models; namespace myop.Controllers { public class Introspect { public bool active { get; set; } public string scope { get; set; } public int exp { get; set; } public int iat { get; set; } public string sub { get; set; } public string aud { get; set; } public string iss { get; set; } } [Route("op/[controller]")] [ApiController] public class IntrospectController : ControllerBase { private readonly myopContext _context; string CLIENT_ID; string CLIENT_SECRET; string TOKEN; public IntrospectController(myopContext context) { _context = context; } // POST: op/introspect [HttpPost] public async Task<ActionResult<Introspect>> doPost() { string body = await new StreamReader(HttpContext.Request.Body).ReadToEndAsync(); string[] p = body.Split('&'); for (int i=0; i<p.Length; i++){ string[] values = p[i].Split('='); switch(values[0]) { case "client_id":CLIENT_ID=values[1];break; case "client_secret":CLIENT_SECRET=values[1];break; case "token":TOKEN=values[1];break; } } var client = await _context.Clients.FindAsync(CLIENT_ID); if (client == null) { return null; } if (client.AccessType == "confidential") { if (client.ClientSecret != CLIENT_SECRET) return null; } else if (client.AccessType == "public") { if (client.GrantType != "password") return null; if (CLIENT_SECRET != null) return null; } else { return null; } string ISS = "http://raspberry.pi:5000"; string SCOPE = null; string SUB = null; string AUD = CLIENT_ID; int IAT = 0; bool ACTIVE = false; var token = await _context.Tokens.FirstOrDefaultAsync(e => e.AccessToken == TOKEN); if (token != null) { SCOPE = token.Scope; SUB = token.Id; int unixTimestamp = (int)(DateTime.UtcNow.Subtract(new DateTime(1970, 1, 1))).TotalSeconds; IAT = (int)(token.Iat.Subtract(new DateTime(1970, 1, 1))).TotalSeconds; if (unixTimestamp - IAT < 60) ACTIVE = true; } Introspect introspect = new Introspect {active = ACTIVE, scope = SCOPE, exp = IAT==0 ? 0 : IAT + 60, iat = IAT, sub = SUB, aud = AUD, iss = ISS}; return introspect; } } }期限切れのTOKEN
$ curl -d "client_id=client1&token=5CCC70961A514EA9A95B195E99B4B754" http://raspberry.pi:5000/op/introspect {"active":false,"scope":"openid","exp":1584998144,"iat":1584998084,"sub":"user01","aud":"client1","iss":"http://raspberry.pi:5000"}期限内のTOKEN
$ curl -d "client_id=client1&grant_type=password&username=user01&password=user01" http:/ /raspberry.pi:5000/op/token {"access_token":"CC8B08A746FC445A98320C7580AB5972","expires_in":60,"token_type":"bearer","scope":"openid"} $ curl -d "client_id=client1&token=CC8B08A746FC445A98320C7580AB5972" http://raspberry.pi:5000/op/intro spect {"active":true,"scope":"openid","exp":1584998352,"iat":1584998292,"sub":"user01","aud":"client1","iss":"http://raspberry.pi:5000"}知らないTOKEN
$ curl -d "client_id=client1&token=CC8B08A746FC445A98320C7580AB592" http://raspberry.pi:5000/op/intros pect {"active":false,"scope":null,"exp":0,"iat":0,"sub":null,"aud":"client1","iss":"http://raspberry.pi:5000"}
- 投稿日：2020-03-24T06:20:09+09:00
dotnet コマンド要約
はじめに
最近dotnetを使い始めたんだけど、どういうふうにコマンド打っていくのか毎回調べている状態で、面倒なのでここにまとめて簡単に記そうと思います。
適当に記してるだけなので、間違っていたら指摘ください。。。基本的なコマンド
プロジェクト開始
dotnet new console -n <プロジェクト名> -o <作成フォルダ名>とりあえず動作確認
dotnet runNuget パッケージ追加
dotnet add package <パッケージ名>尚、基本的なNugetの参照先は
https://www.nuget.org/
のサイトから取ってくる。自己完結型アプリケーション作成
自己完結型：PCにインストールされているFWに依存しないアプリケーション
dotnet publish -c Release -r <RID>RID:https://docs.microsoft.com/ja-jp/dotnet/core/rid-catalog#windows-rids
使用頻度の低いコマンド(所感)
clean
dotnet cleanソリューションファイル作成
dotnet new sln要は、makeみたいなもの。ビルドタスクを回す際は次のコマンドで実行
ソリューションビルド
通常のリリースビルド
dotnet msbuild -property:Configuration=Release自己完結型パッケージのビルド作成
dotnet msbuild -target:Publish -property:RuntimeIdentifiers=<RID>依存関係解決
Nuget等の依存関係の解決コマンド。
ただ、ビルドする時に自動で行ってくれるため基本的にコマンドを明示的に実行する必要性があまりない。dotnet restoreNugetパッケージ作成
Nugetパッケージ作ると他のプロジェクトファイルから参照できるようにNugetパッケージを作成する。（具体的な使い方はよくわからない・・・）
dotnet pack単純なビルド
ビルドも、publishやrunする時に勝手に行われるため明示的にやる必要性は基本的にない。
dotnet build -c Releaseおわりに
とりあえず、なんとなくわかる範囲で調べた内容をまとめました。
なんか追加で調べて行ったら、またここに追加していこうと思います。
- 投稿日：2020-03-24T00:42:00+09:00
C#でもスクリプト実行したい！
はじめに
C#er諸兄は人生の中で一度は、
C#をインタプリタ言語のように扱えたらどれだけよかったろう...と考えたことがあるだろう。
それ、できますよってお話。
導入
前提
.NET Core 3.x以上が導入されていることが望ましい
- もちろん
C#er諸兄は問題ないよなぁ？！
VSCodeが導入済み かつC#プラグイン導入済みdotnet-script のインストール
以下のコマンドを実行する
powershelldotnet tool install -g dotnet-script以上、終了
使い方
その１：シェルで使う
サンプルは
Powershellを利用しているが、dotnetコマンドが使えれば別になんでも良い。
以下のコマンドを実行する
powershelldotnet scriptインタプリタ環境が起動するので、あとはお好みで
powershellConsole.WriteLine("Hello, World!!");実行結果
ただ表示したいだけなら
Console.WriteLine()すらいらないpowershell"Hello, World!!"実行結果
もちろん変数に保存した値の確認も同様に可能
powershellvar msg = "Hello, World!!"; msg実行結果
.NETと言えばNuGet. もちろん対応powershell#r "nuget: Utf8Json" var json = @"{ ""key"":100 }"; var data = Utf8Json.JsonSerializer.Deserialize<dynamic>(json); data["key"]実行結果
その２：
VSCodeで使う / 単一ファイル
適当なところに
.csxファイルを作成して、VSCodeで開くpowershellNew-Item main.csx code ./main.csxとりあえず
Hello, World!!をするmain.csxConsole.WriteLine("Hello, World!!"); // .csxだと直値を出力とかはできない. 以下はエラーとなるので注意. // "Hello, World!!"実行してみる
powershelldotnet script ./main.csx
NuGetも利用可能main.csx#r "nuget: Utf8Json" var json = @"{ ""key"":100 }"; var data = Utf8Json.JsonSerializer.Deserialize<dynamic>(json); Console.WriteLine(data["key"]);実行結果
その３：
VSCodeで使う / インテリセンスを効かせるその２の方法でも十分と言えば十分だけれども、
インテリセンスのない C# なんて C# じゃない！！
ということで、インテリセンスを効かせた方法も紹介しておく。
適当なところにフォルダを作成して、
C#スクリプトプロジェクトを作成するpowershellmkdir ./work dotnet script init code .実行結果
インテリセンスが効く！
NuGetパッケージもVSCodeを開きなおせばインテリセンスが効くように！！main.csx// 以下を追記して、VSCodeで開きなおす #r "nuget: Utf8Json"実行結果
※インテリセンスが効き始めるまでに時間がかかる場合がありますもっと詳しく！
公式サイト を読み込むのだ！！
おわりに
私は
dotnet-scriptは、これでがっつり開発！というよりか、わりと断片的なコードの動作確認のために利用しています。
スクリプト利用するならF#の方がかなり優秀だと思います。ただ、実用に耐えうるレベルのものなのでドンドン使っていってもいいのでは？と思います。
以上、閉廷！！
