光波測量結果とグーグルアース画像とを平面座標系で補正する試み　aim for suitability correction btw electro-optical distance measurements and Google Earth’s images using the coordinate plane

1. Google Earthとのつきあい方

　光波測距儀設置点をGEPro情報から捉える：後方交会法では，後方交会法の限界を思い知った結果になった。ドローン測量の際に，公共基準点を捉えられなかったことが今なお，響いている。現地に再度行った方が，よほど人生を生かすことになるはずだが，だらだらと不完全なドローン測量結果を生かすべく，闘っている？
　ドローン測量と併せて光波測量を実施したが，光波測距儀の設置位置が数メートル規模で確定できない。gpsはこの精度では信頼できない。それではと，Google Earthを使おうとのめり込んだ。この限界については，グーグルアースプロの限界にぼくが得た情報をまとめている。その結果，Google Earth上で何らかの作業をするのは，狭い範囲の測量情報を取り扱うのは，時間をドブに捨てるのと同じだと，遅ればせながら，わかった。GEが提供するプログラムがあるが，これは児童のお遊びツールとしてのみ，何とか価値があるかも知れない。

2. 平面上での作業を選択

　現在の悩みは，前述のように，もともとのドローン測量の周辺的な問題もあるが，Google Earthの画像と光波測距儀のgcpのズレを気にしてしまったことにある。Google Earthの画像が当方の使用目的には不適当であることがわかっても，ドローン測量から得られたortho.tifをGoogleアースに落とした際にほぼGoogle Earthの画像と一致していて，違和感を感じさせないということに意味があると，グーグルアースの画像がモザイク写真であっても，思っている。
　このページ作成の前に，前述の後方交会法のページで最後に掲載したGoogle Earthの画像から後方交会法で得たプッシュピンを削除し画像を図１に示す。赤いプッシュピン（電柱byKなどの個別名称）はGoogle Earth上で，現地調査時の現場写真やドローン測量で得られた空中写真に基づいて，手作業で落としたものである。これをここでは仮にGoogle Earth系A〜Fとする。
　茶色のプッシュピン（a-GE〜g-GE，-GEを付加しているのは光波測距儀の設置位置を現場写真なろからGE上に落としたことに由来する）は光波測距儀で得られたものである。これをここでは仮に光波系a〜fとする。

　今後の作業過程を次に。

ステップ1　光波系a〜f（計7点）とbase-GE（GE画像にプロットした光波測距儀の想定設置位置）（以上，計8点。平面直角座標系1の形）。Google Earth系A〜F（計7点）についてはkml出力して，経緯度（十進角）を平面直角座標系1に一括変換する。

Wikipedia由来か。次の情報を得て，今後は，十進角で表現したい。
「十進角（じゅうしんかく、英：Decimal degrees (DD) ）とは、角度の表現法の一つである。緯度と経度の地理座標を度分秒法が北緯a度b分c.d秒、東経e度f分g.h秒と表現するのに対して、十進角ではa.bcd度、経度e.fgh度といったように分秒の部分を小数、十進法で表す。」

ステップ2　この二つの系を一つの散布図にプロットして画像出力して，それをイラストレーターに配置する。この一連の作業の目的はGoogle Earthのモザイク写真との整合性を取ることなので，光波系a〜f（計7点）のみ，トレースしてベクトル化して，手作業？で回転および平行移動をする。

　なお，この回転と平行移動は現実的でなければならない。平行移動は最大2mほどだろうし，回転は数度以内であろう。光波測距儀では回転軸はコンパスで実行されるが，大山の自衛隊基地のレーダーの影響が考えられるので，数度の回転を要する可能性はある。

ステップ3　ステップ2の結果に基づいて，幾何学的な数値計算を実施する。

ステップ4　その結果に基づいて，光波測量の結果を回転および平行移動を実施する。

ステップ5　ステップ4で得られた結果を十進角に変換してKMLファイルを作成してGoogle Earthに取り込む。

ステップ6　Google Earth上で，ノータッチのGoogle Earth系A〜F（7点）と変換後の光波系a〜fとベース（8点）の整合性を一応，確認する。この作業を通じて，当然，みかけ状の改善は見られるはずである。

3. 計算作業

3.1 ステップ１

　Google Earthから出力したKMLをみると，図2のように，

経緯度値が2セット見られる。
<LookAt></LookAt>内に，128.5300162650217 27.35264129535775
<Point></Point>内に，128.5300076128268,27.35263507066106
KML リファレンスには，<LookAt></LookAt>について，次の説明がある。

Feature から派生する要素に関連付ける仮想カメラを定義します。LookAt 要素は、表示対象のオブジェクトに対する「カメラ」の位置を設定します。Google Earth では、ユーザーが [場所] パネルのアイテムや 3D ビューア内のアイコンをダブルクリックすると、「空を飛ぶように」視点が移動し、LookAt で指定した視点のビューに切り替わります。

というわけで，建物南西角の場所は，<Point></Point>内の経緯度を採用する必要がある。

　国土地理院の平面直角座標系への一括変換では，DMS形式の形に限定されるので，十進角からDMS形式にまずは変換しなければならない。経緯度座標系↔平面直角座標系　にその手法は掲載している。bl2xy.inのファイル内には小数点以下4桁であるが，10桁でも問題なく，計算される。bl2xy.outの平面直角座標系の座標値は小数点以下4桁までに限定される。

表１　Google Earth系A〜F（7点）の経緯度（十進角）から平面直角座標系1へ

3.2 ステップ２

　エクセルでの光波系a〜f（計7点）とbase-GE（GE画像にプロットした光波測距儀の想定設置位置）（以上，計8点。平面直角座標系1の形），と，Google Earth系A〜F（計7点）の散布図を作成する。

　図3にエクセルで作成した両グループの散布図を作成し，左の修正対象の光波系a〜f（計7点）とbase-GE（計8点）の散布図を，右の固定とみなすGoogle Earth系A〜F（計7点）の散布図にコピペした。

　図3右の散布図では，固定側のマーカーを赤い十字マーカーとしている。右の散布図を基図に，イラストレーターでズレの傾向を示す。図3右の散布図を見て，意外とズレが小さいので，方針を修正した。
　エクセルの散布図をイラレにコピペすると文字化けする。MS ゴシック，Osaka，Times New Romanに変更すると回避できたという情報があるが，ぼくは全然だめ。Mac Sonoma 14.1での，Microsoft Excel 16.79.2 から，Adobe Illustrator 2023 へのコピペだ。ひょっとするとと，Adobe Caslon Pro (Regular)にしてみると，見事に文字化けしていない。
　ついでに言えば，以前のマックでは，エクセルからコピーした図表はすべて，自由に分割，削除ができた，そういうマッキントッシュがあった。今は，ユーザーの使い勝手無視の体制だ。アプリ間の連動性が悪い。マイクロソフトとアドビの間でのやり取りはないのだろうか。いずれもマックアプリメーカーとしては老舗なのに。

3.3 ステップ３

　図４の右図はイラストレーターで描いたものである。(base)を頂点として生まれる誤差三角形は細長く，北西部分はF(f)を，南東部分はG(g)の例を描いている。F(f)，A(a)では，グーグルアース画像に比べて，(2.998 + 2.559)/2 = 2.7785º，G(g)，B(b)では，(2.251 + 1.972)/2 = 2.1115º，となる。

　Google Earth系A〜F（計7点）は固定点である。光波系a〜f（計7点）とbase-GE（計8点）は修正が可能である。base-GEを移動すると，光波系a〜f（計7点）も計算式に従って，変わる。iterationの手法を使えば，北西部分と南東部分の各2点の誤差角度の平均値を揃えることは可能である。(base)を最大2mほどだろうか南方に移動すれば一致する筈である。

　このプログラムを作成することはできるが，また時間を費やしてしまう。このウェブページの作成の前に予想した通りになった。このページの作業は，モザイク写真で構成されるグーグルアース画像と合わせる作業である。固定点のGoogle Earth系A〜F（計7点）の座標値は，実はモザイク写真と対応しているものではない。ちょっとややこしいく，ぼく自身も吹き飛ばされそうになるが，この座標値はモザイク写真の歪んだ像にプロットしたものである。だから，結局，モザイク写真に振り回されていることになる。もういいんじゃないか，と思う。回転すら必要無いのではないか，と思ったりしてしまう。

　でも，それでは，なーんで，ここまでグーグルアース画像との対応を取る「努力」をしてきたのか，と思う。とにかく，(base)の南進作業はやめよう。図４で得られた，回転角をどう生かすかを考えよう。北西部分の２ペアの平均値，南東部分の２ペアの平均値，そしてまとめて４ペアの平均値，いずれかを採用しよう。図１と図４右図を睨んで，結局，まとめて４ペアの平均値，を採用するのが無難であろう。2.445ºだ。

　光波系a〜f（計7点）は，反時計回りに2.445º，歪んでいるのだから，時計回りに2.445º回転すれば良いことになる。

　さて，図４右図の光波系a〜f（計7点）の，c, d, e, の三点のズレの傾向を確認したい。青い矢印で感覚的にズレを表現しているが，まあ，トレンドとしては悪く無い。　

3.4 ステップ４

　表２は，光波測距儀のgcp読値を整理したものである。この場の磁気偏角は5º19′W。光波鏡筒軸をコンパスで磁北にゼロセットしてターゲットの方位角を求めて行く。例えば，a点の場合，HR読値: 325º9′50″。磁北は真北に対して5º19′西偏しているので，真北は時計回りに5º19′の位置にある。そこで真北からの水平角（表２の「補正済み真北からの角度」）は，325º9′50″ – 5º19′ = 319º50′50″ となっている。

　3.3スタップ３の結果を，表２に反映させる必要がある。HRは光波測距儀の読値だから触らない。rad列の計算式は，「補正済み真北からの角度」の数値を使っている。磁北の西偏を補正するために，「磁北からの読値」から西偏分差し引いたのであるが，これは「磁北からの読値」からすると反時計回りの処理をしたことになる。図４右図の(base)を中心とするして，光波系a〜f（計7点）をいわば，反時計回り↪️に余分に回転した結果と言えるので，5º19′ – 2.445º = 5.317º – 2.445º = 2º52′ に抑えれば良いということになる。2次元でさえ，回転を直感的に理解するのはぼくは苦手だ。

表３　Google Earth画像に光波測距儀位置を求めて，光波測距データを整理したもの

計算結果が次の表２だ。

　この結果を使って，図３を更新した結果が図５である。回転後の光波系a〜f（計7点）は，Google Earth系A〜F（計7点）へかなり接近した。北東隅の電柱トップも接近している。

　図５下段では，baseからの挟角が回転前からすると回転後は小さくなっている。

図５　回転後の光波系a〜f（計7点）はGoogle Earth系A〜F（計7点）にかなり接近した

3.5 ステップ５

　ステップ４までで所期の目的は達成した。グーグルアースで表示する必要性は全くない。
　とはいえ，この結果に基づいてortho.tifが作成されるので，このオルソ画像上にgcpの正確なプッシュピンを表示する意味はある。2018, 19年の調査情報もグーグルアースのオルソ写真上に掲載したいと思う。
　

　回転後の光波系a〜f（計7点）の平面直角座標系から経緯度座標系への変換をこれから実施したい。出来たものが次のようだ。

　次は，表５の枠囲の数値からKMLファイルを作成することになる。

以上，Dec. 8, 2023記。

　表５の十進角の値を使って，kmlファイルを作成し，グーグルアースに落としたのが次の図６である。良い感じだ。2018, 2019年に川口さんと光波測量した際の光波測距儀の設置点も表示したのが図７である。

　回転後のポイントg-rotは，2018, 2019年に川口さんと光波測量した際の光波測距儀の設置点ls_5に一致するが，図７ではずれている。ortho.tifが手に入れば，2018,2019の座標系のls_5をg-rotに平行移動して，その後に必要ならば回転してortho.tif上に2018,2019の座標系を載せたいと思っている。