1. 課題

　フィールドワークでの測量に問題が生じた。2020年測量のdrone３Dマップと，2018年の光波測距儀を使って得られたネットワークが対応しない。droneマップは光波測距儀で得られたgcpとの対応関係が取られている。そのgcpを使ってdroneマップをコントロールしているので問題はない。この場所には自衛隊のレーダーサイトがあって，コンパスが正常に作動していない危惧があり実際に西偏角度が想定されるものと異なっている。

　光波測距儀の設置点は14カ所で個々の設置点で幾つか測点があって，前視と後視で，位置関係を一応把握したのであるが，手抜きしてGPS計測値を光波測距儀の地球座標値として使用した。これが今抱える問題の元凶だ。そこで，2018年の測量時の測点の現場写真からdrone３Dマップと対応できる測点を利用して，ls_1〜_14の光波測距儀の設置点を求めたいと思う。

2. ls_2の設置点を求める

　ls_2の放射星形区で３Dマップ上での２カ所を現場写真から特定できた。その一つは，図１のls2_20であり，もう一つは図2のls2_11である。後者は光波ネットワークとも一致している。それを図3に示す。

　ls2_20の方は，図２に示すようにかなりズレている。図１の現場写真ではa-rotの下に見えるが，図2では3mほど内陸にずれ込んでいる。

　求めたいのは図１の，赤字で「かわぐち」としたls_2の設置点の３Dマップ上の座標値である。

3. ls_2の3Dマップ上の座標値を求める

3.1 計算法

　当初，3D空間での三角形平面を想定して計算すべきと考えた。しかしながら，熟考するとその必要性が無いことが明らかになった。光波測距儀の構造は台座を水平にして，光波測距儀が装着されている。コンパスで磁北をゼロ度として，台座上を時計回りに回転して，光波測距儀をターゲットに向けて，直線距離を求める。その測距儀の水平角と垂直角を記録する。これだけだ。今，問題なのは３Dマップでのls_2設置点を決定することであり，2D空間に過ぎない。

　連立方程式は簡単にできる。文字表現の簡略化したい。現在の設置点をO’点として，3Dマップに載せうる設置点をO点（Xo, Yo）とする。ls_2放射星形区の２点のうちの１点について，現在の位置をA’点として，3Dマップに載せうる設置点をA点（Xa, Ya）とする。他の１点については，現在の位置をB’点として，3Dマップに載せうる設置点をB点（Xb, Yb）とする。

　光波測量の結果から，次式が成り立つ。赤字表現しているのが未知数である。他の数値はすでに得られた計測値または3Dマップから読み取りできる。
|OA|=|OA’|=√((Xa-Xo)^2+(Ya-Yo)^2)　⋯⋯⋯⋯⋯式(1)
|OB|=|OB’|=√((Xb-Xo)^2+(Yb-Yo)^2)　⋯⋯⋯⋯⋯式(2)
　上記の計算式が成り立つ理由は次のようである。光波測量の結果に基づいて位置ベクトルを３Dマップにプロットするのであるが，ls_2の設置点が決まらないだけであって，その測定結果のうち，水平距離，垂直距離，さらには視準線の間の挟角は変わらない。

　式(1), (2)の連立方程式から，(Xo, Yo)の解が２個現れる。A点を中心として半径|OA|の円と，B点を中心として半径|OB|の円との交点は２点形成されるからである。3Dマップを見ている作者からするといずれが正解かは自明ではある。

　式(1)，(2) は何れも円の方程式である。両円は交わって，その二つの交点を通る直線が成立する。その観点で整理したいと思う。式(2)と(3)はそれぞれ，次のように整理できる。
(Xo– Xa)^2+(Yo– Ya)^2 = |OA|^2　⋯⋯⋯⋯⋯式(3)
(Xo– Xb)^2+(Yo– Yb)^2 = |OB|^2　⋯⋯⋯⋯⋯式(4)

　二つの交点を通る直線の式は，式(3) – 式(4)，で次のように，得ることができる。
2(Xa-Xb)Xo + 2(Ya-Yb)Yo +〔(|OA|^2 – |OB|^2) – (Xa^2 -Xb^2) – (Ya^2 – Yb^2)〕= 0⋯⋯式(5)

　両円の交点の計算は，式(3)と式(4)の何れかの円と式(5)の直線との交点を求めれば良い。ここでは，式(3)の円と式(5)の直線の交点を求めることにする。
　式(3)の，中心座標（Xa, Ya），半径 |OA|の円，と，式(5)の直線の交点をM, Nとして，交点を表現したいのであるが，かなり複雑になるので，式(5)を簡潔化して次に表現する。

e=2(Xa-Xb), f=2(Ya-Yb), g=〔(|OA|^2 – |OB|^2) – (Xa^2 -Xb^2) – (Ya^2 – Yb^2)〕として式(5)は次のようにする。
eXo + fYo + g = 0⋯⋯式(6)

以上，ここまで，Jan. 25, ’24記。

　さて，Microsoft Excelの計算式の書式で，交点MNの座標値を次に示す。

Xm=(e*D–f*sqrt((e^2+f^2)*|OA|^2-D^2))/(e^2+f^2)+Xa⋯⋯式(7)
Ym=(f*D+e*sqrt((e^2+f^2)*|OA|^2-D^2))/(e^2+f^2)+Ya⋯⋯式(8)

Xn=(e*D+f*sqrt((e^2+f^2)*|OA|^2-D^2))/(e^2+f^2)+Xa⋯⋯式(9)
Yn=(f*D–e*sqrt((e^2+f^2)*|OA|^2-D^2))/(e^2+f^2)+Ya⋯⋯式(10)

　なお，D=abs(e*Xa+f*Ya+g)⋯⋯式(11) とする。

　式(7)と式(9)，そして，式(8)と式(10)，の違いは第二項と第三項の間の符合に見える。

この計算プロセスを踏まえたMicrosoft Excelファイルは次節に用意する。

3.2 測量値などの準備

　ls_2放射星形区で，２カ所 ls2_11と_20の3Dマップ上の位置決めができている。ls2_11の座標値は変わらず，ls2_20だけを現場写真（図１）と対照して決定する必要がある。図4で特定している。光波測距儀による測量結果では，altitudeは1.171mだから，図4の右上に示した表では二行目が近い。まあ，これよりも海側下方であろうが，ほぼ対応しており，図4右上表の，Index 107222891を採用したい。
X: -96035.546875, Y: -625538.500000, Z: 1.243085

　以上から，式(3), (4)の既知数値は次のようになる。なお，注意すべきなのは世界測地系の (x, y)=(northing, easting)と3Dマップの(X,Y)=(easting, northing)は入れ替わっているということである。学校教育で使用するデカルト座標系に替わっている。繰り返すが，ls2_11の座標値は3Dマップと一致している。両円の半径は光波測距儀の読値HD horizontal distanceに対応している。

ls2_11を中心とする円関係：　(Xa, Ya)=(-96021.025, -625506.207), |OA|=42.830⋯⋯式(12)
ls2_20を中心とする円関係：　(Xb, Yb)=(-96035.547, -625538.500), |OB|=65.508⋯⋯式(13)

式(12), 式(13)は式ではないが，便宜上，式番号を与えた。「3.1 計算法」の計算プロセスをMicrosoft Excelに埋め込もうと思う。　

3.3 Microsoft Excel計算シート

　計算シートを作成し，実行した。しかしながら，ls2_20の3Dマップ上の認定を失敗したようで，求める結果を得ることはできなかった。そのプロセスは次の節にしめす。Microsoft Excelの計算シートでは，まずは次のテキストボックスを配置した。

————————————————

解　説

　３Dマップ上で特定された２地点から光波測距儀の設置点を求める手法をここに示す。

現在，設置点としてきたものをO’点として，3Dマップに適切に載せうる設置点をO点（Xo, Yo）とする。この例では，ls_2放射星形区の２点のうちの１点について，現在の位置をA’点として，3Dマップに載せうる設置点をA点（Xa, Ya）とする。他の１点については，現在の位置をB’点として，3Dマップに載せうる設置点をB点（Xb, Yb）とする。

　結局のところ，A点（Xa, Ya）を中心とする半径|OA|の円と，B点（Xb, Yb）を中心とする半径|OB|の円の交点を求める問題になる。連立方程式は次のようになる。

(Xo- Xa)^2+(Yo- Ya)^2 = |OA|^2　⋯⋯⋯⋯⋯式(1)
(Xo- Xb)^2+(Yo- Yb)^2 = |OB|^2　⋯⋯⋯⋯⋯式(2)

　交点MNの座標値を次に示す。

e=2(Xa-Xb), f=2(Ya-Yb), g=〔(|OA|^2 – |OB|^2) – (Xa^2 -Xb^2) – (Ya^2 – Yb^2)〕⋯⋯式(3) として，

Xm=(e*D–f*sqrt((e^2+f^2)*|OA|^2-D^2))/(e^2+f^2)+Xa⋯⋯式(4)
Ym=(f*D+e*sqrt((e^2+f^2)*|OA|^2-D^2))/(e^2+f^2)+Ya⋯⋯式(5)

Xn=(e*D+f*sqrt((e^2+f^2)*|OA|^2-D^2))/(e^2+f^2)+Xa⋯⋯式(6)
Yn=(f*D–e*sqrt((e^2+f^2)*|OA|^2-D^2))/(e^2+f^2)+Ya⋯⋯式(7)

なお，D=abs(e*Xa+f*Ya+g)⋯⋯式(8) とする。————————————————

　図５に，この計算部分を掲載する。

　図５のダウンロードファイルは用意しているが，ls_2放射星形区内の２点の選定について，試行錯誤の状態のため，成功の後のファイルを掲載したい。

3.4 計算結果の検証１

　図５の計算結果は，Ⅳ〜Ⅶに見えるが，位置関係から見ると，(Xm, Ym)を採用できる。内積計算でも評価できるが省略する。当方のMicrosoft Excelファイル Okinoerabujima川口木庭_2018Jun14 参考ファイルJan2024.xlsxの沖永良部島測量ls No.2では，ls_2の光波設置点はGeographicaで得られたGPS表示を採用しているが，その数値を，上記(Xm, Ym)に替えた。このスプレッドシートは，northing X, easting Yに設定しているので，図5の座標定義と逆になっているので注意が必要である。

　 ls_2の設置点を上記(Xm, Ym)に替えたのあるが，適切な結果が得られなかった。光波の視準線の磁北設定も変更する必要があった。もちろん，この視準角のゼロ位置を変更する必要性は当初から考えていたが，結局，3Dマップの真北設定でのXY座標値の土俵に入っているのだから，光波の測定値も磁北の偏角を考慮する部分は削除する必要があったのである。

　「図６　光波計測値の表に図５の計算結果を入力して磁気偏角をゼロに」を参照頂きたい。この計算過程を実施する前に，もともとのls_2の計算シート（沖永良部島測量 ls No. 2）をコピーして，この右隣に配置した。そして，最上行の「ls-site northing X」, 「ls-siti easting Y」カラムに，上記(Ym, Xm)を入力し，さらに，磁気偏角をゼロにするために，HRの角度を，補正済み真北からの角度，にそのまま，コピペした。その変更を赤字で示している。Loc. 11はls2_11に，Loc. 20はls2_20に該当している。もともとのls_2の計算シート（沖永良部島測量 ls No. 2）の値と比べると，ls2_11は，(x North, y East) = (0.19, 0.215)，ls2_20は，(x North, y East) = (0.46, 1.901)，となった。

　この結果から，次の表に示したcsvファイルを作成し，CloudCompareに読み込んだ。

site_name	X easting	Y northing	altitude	Red	Green	Blue
ls2_5	-96002.656 	-625419.660 	1.872 	0	0	255
ls2_6	-96018.265 	-625483.024 	0.466 	0	0	255
ls2_7	-96022.656 	-625483.303 	0.153 	0	0	255
ls2_8	-96024.513 	-625483.271 	-0.493 	0	0	255
ls2_9	-96032.904 	-625480.370 	-0.113 	0	0	255
ls2_10	-96015.073 	-625486.809 	2.183 	0	0	255
ls2_11	-96020.810 	-625506.010 	2.159 	0	0	255
ls2_12	-96040.427 	-625525.617 	-0.510 	0	0	255
ls2_13	-96039.064 	-625525.209 	0.437 	0	0	255
ls2_14	-96032.476 	-625524.511 	0.432 	0	0	255
ls2_15	-96030.971 	-625523.683 	0.917 	0	0	255
ls2_16	-96029.524 	-625529.911 	2.325 	0	0	255
ls2_17	-96039.251 	-625539.972 	-0.491 	0	0	255
ls2_18	-96037.813 	-625540.109 	0.665 	0	0	255
ls2_19	-96034.701 	-625539.333 	0.733 	0	0	255
ls2_20	-96033.646 	-625538.043 	1.171 	0	0	255
ls2_21	-96033.215 	-625538.873 	2.464 	0	0	255
ls2_22	-96032.121 	-625554.452 	0.614 	0	0	255
ls2_23	-96026.456 	-625577.528 	0.404 	0	0	255

　図7を見ると，ls2_20などの位置はより左手（つまり海側）に移動しているが，例えばls2_18は，ls2_17付近に位置する必要がある。図8は新プロットだけを表示したものであるが，垂直関係に問題はないようであるが，図7のls2_19は本来，ls2_18が載っている緑藻の面の陸端に配置すべきものであり，光波測点が3Dマップとうまくフィットしていないことが見えるのである。

　結果からすると，ls2_20の3Dマップでの特定に誤りがあったようである。現場写真を見直して，ls2_19の方が特定しやすいように思われるので，これまでls2_20で実行したことを，ls2_19で試みることにした。

以上，14:36，Jan. 28, 2024記。

3.5 計算結果の検証２

　ls2_19の図10での想定位置位置を，図９の現場写真から決定した。図10の表の2の位置が適当と思われる。図10の鮮明度は低く，2階調にしている。このデータは次のようである。
X easting =-96036.617188, Y northing = -625539.000000, Z = 0.711838

　ls2_11をCloudCompareで観察する限り，大きなズレはみにくいが，なにせ，drone撮影の影になっており，特定するのはかなり難しい。ls_19を採用してもうまく行かない場合は，ls2_11の代わりをみつける必要性がある。大変だなあ。

3.6 計算結果の検証3

　さて，これまで触らなかったls2_11に代わって，ls2_6（Z=0.466 m a.s.l.）の現場写真を観察して，3Dマップ上の座標値を求めた。図12右上のpicked points list表には２点が見えるが，下方の方を採用した。
　なお，図12のls2_6の海抜高度は高いように見えるが，Cloudを透過して表示されているからである。list picking作業で，水平位置が海側に移動することになる。

X easting = -96021.109375, Y northing = -625482.125000, Z = 0.445537である。

　このページは一昨日以来，触らなかった。ls2_6付近の表示の差を図13，14に示している。図2では北縁の測点はls2_11までしか見えていない。図13, 14の場所は図2よりも北方にある。図11のプリズム三脚の位置がls2_6にあたるが，図13では3Dマップ上のls2_6ポイントは，図11よりも更に内陸側に位置してる。図14を見ると，ls2_6は，いい按配に見える。それゆえ，両円の中心として，半径|OA|を作るA点をls2_11からls2_6に替えることで，大幅に改善されたことがわかる。

　このように，A点=ls2_6，B点=ls2_19,とすることで，図14のように，かなり改善を見た。では，ls2_19付近はどうなったのか。図7と比較すると，図15ではより改善され，光波の測点はより海側に移動した。

　A点をls2_11からls2_6に替えることで，ls2_6付近も，ls2_19付近も改善されたのである。ぼくの性格から，ls2_20とls2_6の組み合わせはどうかとやってみたが，劣化した。今のところ，ls2_6, ls2_19の組み合わせがいいようだ。

　図15で気になるのは海側にずれ込み過ぎていること，図16で気になるのは，ls2_6. _7, _8が左方にずれ込んでいることである。

　ls2_19付近の3Dマップと現場写真の対照はたぶん，これ以上の改善は難しい。ls2_6付近はどうか。これは悩んだ上での決定であったが，droneで撮影した空中写真の広域の方で見えるのではないか，と考えている。明日は，空中写真を判読して3Dマップと対照して，ls2_6の位置決めをしたいと思っている。

以上，23:18，Jan. 29, 2024記。

3.7 計算結果の検証4

　現場写真から，ls2_6をほぼ特定した。この海側（写真の上方）に波のスプラッシュが見えるが，これは図16のls2_7（Sw-spP）に対応するだろう。ls2_6は，Lw-spPだ。宇和美崎狭域3Dマップで，この場所を特定する必要がある。宇和美崎広域写真で特定したls2_6の位置だ。spratidal platformの陸端だ。

　図18には，宇和美崎狭域3Dマップ上で新たに特定したls2_6は，Point #1にあたる。

X= -96021.5703125，Y= -625484.5625, Z= 0.445600688457

　やっと得たls2_6が新たなA点となる。B点（ls2_19）はすでに確定している。

　両円の交点をもとめ，それをls_2の設置点としてls2_5〜_23の（X, Y）を再計算してcsvファイル ls_2star_ls2_new6_19.csvを用意して，CloudCompareに読み込んだ。図19では，ls2_6付近のもので，右寄りのものが新たな結果で，フィールド写真とよく合っている。

　図20では，ls2_19付近のプロットで，過去で一番改善されたのであるが，いまだ，多少，海側にずれ込んでいる。

　これで，ls_2放射星域については終了しても良いが，ls2_19の位置決めを空中写真判読でも確認したいと思う。

以上，Jan. 30, 2024記。

3.8 計算結果の検証5

　ls2_19を写真（図9）を中心に，宇和美崎写真広域355, 356の実体視を加えて，検証した結果，あらたな位置を特定できた。図21には，新たにls2_19が特定されている（図中のPoint #0）。これまでの最新のプロットと対照すると，ls2_19などがより内陸側に移動できそうである。X= -96036.671875, Y= -625539.3125, Z= 0.754368662834となる。光波でのZ値は0.733m a.s.l.である。

１　ls2_new6とともに，_new19が更新され，両円の交点を新たに得ることができる。
２　得られた交点をls_2光波設置点としてls2_5〜_23の（X, Y, Z）を求める。
３　Microsoft Excelシートからcsvファイルを作成。
４　CloudCompareで読み込んで宇和美崎狭域3Dマップ上にプロットする。

　この結果を次に。図22には，平行に配置されている同名のものの右下が新たなls2_19によるもので，本来，ls2_19であるべき場所のPoint #0には到達できていない。何故だろう。

　 図23の結果は，図21の結果と極めて類似している。平行にズレている。さて？

以上，Jan. 31, 2024記。

4. 計算法を確認したい

4.1 問題の所在

　これまで述べてきたように，図22のPoint#0は，現場写真，空中写真を見つめて，3Dマップにls_19の位置と決定した場所に，CloudCompare上でPicking points list機能を使って，プロットした地点である。ベクトルOBの関連では，Oがls_2設置点で，BがこのPoint#0にあたる。両円の交点がOだから，|OB|=66.509mにあたる。光波測量の計算表に示しているように，|OB|は，HD horizontal distanceに当たっている。そういうように設定した訳だから，交点とls2_19の距離は，|OB|=66.509mの筈であり，宇和美崎狭域3Dマップ上のPoint#0と光波計算表のls2_19は一致する筈である。なのに，なぜ一致しないのだろうか。

　改めて，宇和美崎狭域3Dマップ上に描いた，ls_2光波設置点（図24のls_2setting）と，ls2_5〜23の分布を見たい。

以上，Jan. 31, 2024記。

　両円の交点を求めるこれまでのプロセスを，図25のように図化した。何故に今採用している計算プロセスを経て得られたls2_19が現場写真から特定したPoint #1と一致しないのか，それを特定するためである。

　これからのプロセスを確認すべく，流れを次に示す。

確認１　図25での両円の半径が正しく光波測距儀のHDと一致しているか。
　図25を評価する際に二つの面から調べる必要がある。点A, Bの座標値は，現場写真と空中写真から特定して3Dマップにpicking pointsの座標値にあたっている。これをここで，再掲したい。そしてMicrosoft Excelで計算して|OA|，|OB|が光波のHDに一致することを確認する。
　両円の交点を求めてそれをls_2設置点の座標値として，光波の計算プロセスで得られたls2_6, _19の座標値を使った場合，|OA|，|OB|が光波のHDに一致するか。まあ，一致しないだろう。しかし，一応，この結果を確認したいと思う。

確認２　１の検証結果は，光波の計算プロセスで歪みが生じたということになるだろう。では原因は何か。

確認３　まだやってないからわからないけど，まずは3Dマップが真北を捉えていない。gcpを使って適当に配置している。そして光波で得たネットワークを構成する放射星形区も真北がわからない。両方に原因がある。しかし，3Dマップに光波放射星形区を載せようということだから，3Dマップの擬似真北に光波の放射星形区も合わせないといけない。じゃあ，どうするか。
　図25のベクトルOAの時計回りの角度は，点ls_2settingの座標値を（Xo, Yo）として，現場写真と空中写真から特定された点ls2_19photoの座標値を（Xb, Yb）として，直角三角形の，角 [（Xb, Yo）（この頂点が直角）-（Xo, Yo）-（Xb, Yb）]を余弦定理から求めれば良い。その角度と光波測量時のHR（時計回りの回転角）を対照すれば，西偏用の角度補正のセルを変更すれば，ls_2設置点からの測点の(X, Y)座標値を得ることができるはずだ。

　楽しみだ。今日はタニハで疲れた。明日にした方がいいな。

以上，0:20，Feb. 2, 2024記。

　次にMicrosoft Excelで検証した結果を示す。
————————————————
1. OAとOBを半径とする円の中心と１交点
X Y Z
A相当：ls2_6photo -96021.570 -625484.563 0.446
B相当：ls2_19photo -96036.672 -625539.313 0.754
光波設置点：ls_2setting -96052.835 -625474.797 -0.235

2. |OA|と|OB|の距離と光波HDとの関係 計算結果 光波HD 評価
sqrt((Xa-Xo)^2+(Ya-Yo)^2) 32.754 32.754 当然，一致
sqrt((Xb-Xo)^2+(Yb-Yo)^2) 66.509 66.509 当然，一致
評価：　両円から交点を求めるプロセスに誤りはない。

3. 光波計算過程を得て得られたls2_new6, ls2_new19
X Y Z
ls2_new6 -96021.186 -625483.235 0.466
ls2_new19 -96037.622 -625539.543 0.733
光波設置点：ls_2setting -96052.835 -625474.797 -0.235

4. |OA’|と|OB’|の距離と光波HDとの関係 計算結果 光波HD 評価
sqrt((Xa-Xo)^2+(Ya-Yo)^2) 32.754 32.754 そうか，やっぱり一致
sqrt((Xb-Xo)^2+(Yb-Yo)^2) 66.509 66.509 そうか，やっぱり一致
————————————————

　両円の交点の計算過程でも，これを使った光波放射星形計算過程でも，|OA|, |OB|いずれも変質していないことがわかった。すぐ上の確認１〜確認３での思考プロセスに問題がないことがわかった。

4.2 3DマップのXY座標系に光波放射星形区のX’Y’座標系を合わせる

　用語が定まらないが，このタイトルで考えたいと思う。

　現場写真とdroneで撮影した空中写真を読み込んで，ls2_6とls2_19を中心とする円の交点を使ってきたのだが，点群の分解能と現場写真の質の点から，ls2_19の信頼性が高く，これを優先して計算過程を進めることになる。

　図25の3Dマップの点群を非表示にして，Yo軸からls2_19軸までの時計回りの角度を求めたい。黄色で表示している△BOIに注目して，余弦定理で計算した。

————————————————
5. 余弦定理で直角三角形BOIの∠Oを求める
辺b=辺OI=|Xb-Xo| 16.163
辺i=辺BO=|OB| 66.509
辺o=辺IB=|Yo-Yb| 64.515
cosO: (b^2+i^2-o^2)/2*b*i 0.243
arccos(cosO) 1.325 radian
arccos(cosO)+p()/2 (Y軸からの時計回り) 2.896 ： 3Dマップ上での時計回りのOB角 radian
ls2_19のHRは166度46分40秒を度表示で 166.778
radianに変換 2.911 ：光波計算表でのHR radian
光波計算表でのHR補正値radian -0.015 radians
radiansからdegreesに変換 -0.841 degrees
degree + minute/60 + second/3600 -0°50′27.6″ ：この値を磁気偏角にして実行する。
————————————————
　上記の計算過程で，光波計算表で，実行した。
————————————————
6. 座標回転を考慮した光波計算過程で得られた地点Bと地点Aの（X, Y）座標の結果
現場写真3Dmap 光波座標回転後 差
OB: ls2_new19 Xb -96036.672 -96036.674 -0.002
OB: ls2_new19 Yb -625539.313 -625539.313 -0.000
結果評価： OB軸で座標回転したから当然一致する
OA: ls2_new6 Xa -96021.570 -96021.066 0.505
OA: ls2_new6 Ya -625484.563 -625482.769 1.793
結果評価： 次の項で述べる。
————————————————

この結果を使って，ls2_5〜_23の点情報のテキストファイルを作成して，CloudCompareに読み込んだ。その結果を次に示す。図27はls2_19付近，図28はls2_6付近である。解像度が低くて見にくいが，いずれも海岸地形の断面図作成用に設置した測点に完全に一致している。

　このウェブページ作成の目的を完璧に達成できているのである。

　ls2_19はベクトルOBを重視して，ls2_6の方も見事に達成できたのである。

　で，やっぱり，気になるのは，ls2_6つまり，OAの方である。計算してみた。
————————————————
7. OA軸についても座標回転の偏差を求めるかどうか
　一つの放射星形区で，現場写真から3Dマップとの関係での２地点を特定して，ベクトルOA, OB を中心とする半径HDa, HDbの円の交点を求めてそれを光波測距儀の設置点とするのであるが，地点A, Bを 特定して初めてこの作業が可能である。ここで示している例２．ではOBから光波測量結果をまとめた のであるが， OAについてもこの作業をする意味はない。地点A,Bでは，Aの方が信頼性が高く，Aの結果から得た偏角を採用して光波測量値を処理するしかない。
　光波のOA, OBの視軸がなす角度は，HRb’ – HRa’から得ることができる。radianでは，2.896-1.817=1.079である。一応，ここでのHRb – HRa を求めてみよう。直角三角形 BOAの∠BOAを次に求めてみる。

地点A, Bの座標値を再掲する。
X Y
A相当：ls2_6photo -96021.570 -625484.563
B相当：ls2_19photo -96036.672 -625539.313

辺b=辺OA 32.754
辺a=|OB| 66.509
辺o=辺AB 56.795
十進角（º）
cosO: (b^2+a^2-o^2)/2*b*a 0.521
arccos(cosO) = (HRb-Hra) 1.023 58.590
挟角の差　(HRb-Hra) – (HRb’ -Hra’) radian -0.056 -3.233
なお，(HRb’ – Hra’)は， 1.079 61.822
OA軸のY軸からの時計回り角度 radian 1.873 107.339
ls2_6の時計回り角度 radian 1.817 104.106

光波計算表でのHR補正値radian 0.056 3.233
：OAを採用すると補正角
————————————————

　以上で，一応納得の行く結果となった。
以上，Feb. 3, 2024記。

5. ls_2放射星形区の結果を使って次はls_1区

　これまでの努力で，ls_2は完成した。これに隣接する，ls_1とls_3にアタックすることになるが，まずは_1を。ただし，今のところ作成できた3Dマップは宇和美崎狭域のもので，ls_1の範囲の全域がカバーされていない。広域の3Dマップは未完成だ。
　放射星形区の間では，必ずしも２点での前視後視を実施していないが，ls_1とls_2の間では実行しており，Locs. 5, 6がそれに該当する。3Dマップとの対応関係はls_2で成功したので，ls_1の光波設置点は，ls2_5, 6の(X,Y)値を使って実行できるはずである。

5.1 新たなls_1の座標値を求める

　前掲の「3.ls_2の3Dマップ上の座標値を求める」，のプロセスをまずは実行することになる。(X, Y)値はls_2で最終的に得られたLocs 5, 6の値をls_1のそれとして代用することになる。HDについては，ls_1からのHDを使用する 。
easting 　 northing HD (ls_1)
ls2_5（ls1_5) 点A -96006.388 -625419.183，72.459 （過ってls_2のもの），72.465 （正しくls_1のそれ）
ls2_6 (ls1_6) 点B -96021.066 -625482.769 ，32.754（過ってls_2のもの），33.157（正しくls_1のそれ）

「3Dマップ上の2点から光波測距儀設置点を求める.xlsx」には問題があった。新しく2円の交点を内積で求める（２円交点EnToEn_revised.xlsx）を作成して両円の交点は２点得られた。Google Earthにプロットした旧来のls1_5, _6，GPSのls_1の配置からみて，Xm=-96053.17945, Ym= -625474.51596をls_1設置点と決定することができた。

5.2 新たなls_1の光軸回転角を求める

　「Okinoerabujima川口木庭_2018Jun14参考ファイルJan2024.xlsx」に，上記で得られた新たなls_1設置点の座標値を，入力した。すでにこのページで議論してきたように，さらに，光波測量の放射星形区の測量結果の座標回転が必要になる。「4.4 3DマップのXY座標系に光波放射星形区のX’Y’座標系を合わせる」，で述べた作業を実施する。△BOAを点O（新ls_1）の回りに回転するが，OA軸，OB軸の何れを使用するかで多少異なるが，現場写真を見ると，未完の宇和美崎広域３Dマップがたとえ後に完成しても，点A (ls1_5)を3Dマップにプロットするのは難しい。そこで，OB軸の回転角度を求めることにした。結果として，△の頂点と辺の名称はほぼそのまま使用できる。使用Microsoft Excelファイルは，「3Dmap上で光波星形区をフィット両円交点と座標回転.xlsx」（前掲のファイル名をより一般化してみた）にあたる。その計算シートを図29に示す。

　第29行（赤線の１）にはには余弦定理を使って求めた十進度104.413が見える。第31行（赤線１’）では，それをDMS表示に変換している。そして，その値をls_1星形区のLoc. 6 (ls1_6) のHR(DMS)で差し引いた結果が，赤線２に示す-02:53:53.6である。DMSの引き算の結果は➕️である必要があり，それを，第21〜27行右手に見えるテキスト枠に示している。第34行はその演算行にあたる。

　この「偏差入力値」を使って，ls_1星形区のLocs. 1〜6の再計算を実施する。

　この後の計算の流れは，Microsoft Excelで度分秒表示角の差を求める ，「3. 光波測量結果に」に示している。その図５を，ここでは図30として再掲した。

　これまで述べてきた手法で，ls_3以降の星形区の光波結果を再計算したいと思っている。

　次の図31はCloudCompareで表示したところである。ls1_1〜_6とls2_5〜_13である。

以上，Feb. 19, 2024記。

evil Google

Google browsersには耐えられない。ここ一〜二年ほどか，種々のツールの利用を促されて，Google Chromeの利用頻度は増してきた。Google Driveは100GBの契約をしており，このサイトのバックアップを取っている。Google ChromeをデフォルトWebブラウザーにも設定してきた。

この一カ月ほどか，Amazonなど種々のサイトにアクセスする際に，ぼくがそれぞれのサイトで設定したパスワードが拒否され，Googleに代用するよう強要される。macのメーラーでgmailを受け取ることができるが，返信ができなくなった。gmailユーザーはそんなこと，お構いなしだ。ぼくのgmailアドレスを教えて，gmailユーザーとは交信している。面倒なことだ。

企業のこの種の機能利用の強制策謀は，独占禁止法違反じゃないのか。日本政府や企業のデジタル技術は，世界でも例を見ないほど低劣だ。なんら，対処していない。官僚たちもこの現状を実体験している筈だ，なのに，なにもしない。というか，できない，能力がない。

国民背番号制が進む。日本のマスコミのレベルも低く，この制度が日本の人材レベルでは，外部からの侵入に対処できない。河野デジタル相もえらそうだけど，全然，想像力も能力もない。腐りきった官僚の言うままだ。

最近，イギリスの郵便局システムの問題が再燃している。その元凶は日本の富士通のシステムだ。富士通は，全く能力がない。入社時にはそれなりの人材が入るが，富士通の企業風土が行政べったりで，安楽に生き延びてきたために，社員が努力しない。社内の流れに載っているだけだ。こういうぼくの思いはこれまで滞在してきた大学での富士通システムの低さに呆れてきた歴史がある。流石に官僚が富士通に名指しで，お勉強をせよ，と迫っている。致命的な実害があるからだ。

ブラウザーに戻るが，Googleはひどい。これから脱するのは，ヨーロッパが安いロシア産の天然ガスから離脱する努力に匹敵するかも知れない。アップルもモノポリー企業で，自らの製品をヴィンテージだと称して，サポート体制をサッサと引き上げる。ぼくが所有する古いmacではsafariは使えない。Firefoxでは，最新のテクノロジーが何の問題もなく使えている。

改めて，GoogleとFirefoxを確認もあって比較してみたが，Firefoxには強制力は全く無い。独占企業のドブ臭さは感じない。ドーネーションをしないとなあ。数回はしてると思うけど，本気で。

以上，Jan. 21, 2024記。

さくらのレンタルサーバー

　今日，Jan. 31, 2024に，この問題への対処法が届いた。凄く面倒でDKIM署名に秘密鍵が必要で，この対処法では，使う気にならない。次に，さくらの，解説のリンクを。
DKIM署名、DMARCを設定したい

以上，Jan. 31, 2024記。

1. 自然環境を無視した住宅開発

　日本では，とくに，明治維新以降，西洋の高い土木技術の導入と人口圧や都市化の要請のなか，本来住めない場所でも，住宅開発が積極的に進められてきた。その一つが，犀川開析谷である。金沢市のWebサイトを参照したが，コンテンツの質は他の行政区同様であった。ハザードマップの質は日本のどこでも低い。

2. ハザードマップから

　さて，金沢市のハザードマップを参照して，過去の浸水実績を見たのが次の図１である。図１右手はご自宅の場所をGoogle mapから特定して図１のハザードマップでの位置を同定している。浸水域がモザイク状になっているのは，報告書を参照していないが，排水溝のトラブルによるものであろう。つまり，破堤によるものではない。犀川の水路幅が小さくて上流の域の排水を運ぶような機能が期待されていない。犀川は，川の名前が付いているが，自然の川ではなく，単なる水路になっている。では，安心していいのだろうか。

　図２は，ご自宅の想定水害水位である。0.5m以上〜3.0m未満，河岸浸食となっている。さて，図１とどう違うのか。図１は比較的最近の実績である。リスクは図２をみた方がいいだろう。破堤の可能性は当然ある。その破堤をした際のリスクを図２は示している。すぐ家のそばの河岸が破壊されるということではなく，上流，場合によっては，下流であっても，破堤した場合，0.5以上，3.0m未満の浸水があるということだ。

　それゆえ，集中豪雨などがあって，非難注意報が出た場合，逃げた方が良いことになる。津波について金沢市のハザードマップは，かなり海岸付近に限定されているような図しかないが，10mぐらいの津波が来たら，この水路を遡って，超波することはありうる。ただ，10mの津波が生じるかどうか。5mぐらいの津波だったら，ここまで上がってくるのは難しいかも知れない。

　なお，逃げる場所だけど，行政（業者）の指定先は，全く間違っている。信用してはいけない。今回の能登半島の津波の際の指定避難所に到達できず，津波に呑まれている。東日本大震災の場合も避難場所は，海岸線のそばの海岸砂丘上であった。東北大の津波研究者の推薦場所だった。責任取らずにのうのうとTVにも今なお出ている。ハザードマップや行政の避難場所を信じてはいけない。

3. 破堤の場合のお宅に係わって想定されるルート

以上，Jan. 16, 2024記。途中です。

1. このページ作成の意味

　もともとの関心はGISソフトでの３Dクラウドとgcpの重ね表示の可能性を探ることであった。ここに至る経緯は次のページに示している。

QgisとGrassGISの何れが3Dメッシュと点群ファイルを生かすことができるか　which is better Qgis or GrassGIS for the display and analization of 3D mesh and cloud data ?

　このページの作成過程で，gisソフトでは実現しないことがわかった。MeshLabに期待したが，これも重い3Dファイルを受け付けない点で問題があった。結局，CloudCompareに頼らざるを得なかった。CloudCompareではMeshLabでどうにもならないものでも何とか作業ができるのである。

　で，問題はレイヤー構造で使用出来るかであるが，Danielの情報で，使えることがわかったのである。ただ，図１のように，Plugins: Cloud layers，は灰色になっている。

　Help > About Plugins，で開いた画面が図２である。Cloud layersを選択するが，これ以上進めない。Referencesとして，サードパーティのWiggins Tech website　が見えるが，ここをクリックするとこの注記が消えてしまう。

2. Plugins: Cloud layersのインストールに向けて

　このサードパーティのリンクをCloudCompareを使用しているWindowsで立ち上げた。ドイツ語が表示される。言語を換える選択肢があり，日本語に替えた。その画面が図３のものである。このサイトは2016年11月から始まっている。図３の「Cloud Layers」プラグインをクリックするとmp4のムーヴィーが始まるがCCもなく，音声も出ない。

　図３の表示のすぐ下には，ダウンロードリンクがある。CloudCompareのミラーサイトの認証を受けているようだが，ヴァージョンが2.12.0のままだ。更新されていない。CloudCompareはフリーだけどサポートが要るというメッセージが色つきになっている。

　次の情報がある。ただし，当方の　Editメニューでは，’Edit > Cloud > Create single point cloud‘しか有効でない。Edit > Scalar fields についても下に示された選択肢はすべて灰色になっているのである。

Here are the new menu entries:

• ‘Edit > Cloud > Create single point cloud‘

  • to create a cloud with a single point (set by the user)

• ‘Edit > Cloud > Paste from clipboard‘ (shortcut: CTRL+P)

  • to create a cloud from ASCII/text data stored in the clipboard

• ‘Edit > Scalar fields > Split clouds (integer values)‘

  • will split the cloud in multiple sub-clouds based on the (integer) values of its active scalar field.
  • to be used with a classification SF typically.

• ‘Edit > Scalar fields > Add classification SF‘

  • shortcut to the ‘Edit > Scalar fields > Add constant SF‘ method.
  • creates a ‘Classification’ SF with a constant (integer) value.

　このWebサイトに入ると，ぼくはすごく怪しい印象を持った。どこにもプラグインのダウンロードのリンクはなく，ただ，コンタクトせよという。しかも種々のページを見ていると，bluetoothに繋げようとする。

　Danielさんが，フォーラムで，
Maybe you can use the new ‘qCloudLayers’ plugin?
Another way, maybe more simpler, is to use the ‘Edit > Scalar fields > Split cloud (integer values)’ method. This will create as many (sub) clouds as classes present in your cloud.

という回答をしている理由もそこにありそうな気がしている。

3. 最新バージョンをインストールしてみる

　信頼できるGitHubギットハブで確認した。Version 2.12.0には，確かに，上記の魅力的なプラグインが新たに装備されている。では，何故，最新のステーブル版では灰色になってしまったのか。これがまずは最初のポイントだ。ver. 2.12.0のコードネームKyivを獲得したのもキーウの貢献ゆえだろう。

　V. 2.12.1〜.4まで，特に当該プラグインがinactiveになったのかは，記されていない。

　ではと，最新のアルファ版をダウンロードしてみるのがいいのだろう。GitHubよりももちろん，本家からがいいだろう。本家をみると，より進んだβ版になっている。
CloudCompare　3D point cloud and mesh processing software　Open Source Project

　インストーラーは265MB以上ありダウンロード継続の表示が見えないので非常に取り扱い難い。圧縮版は64MBで軽く，7zで解凍して，ドラッグアンドドロップで，Program Filesに投げ込んだだけだ。ショートカットを作って，デスクトップに配置した。macの方は不思議なのだがダウンロードが始まったのは見えたが，時間をかけて見ると，消失している。２回繰り返したが同様の現象が見られた。ファインダーで検索しても見当たらない。図４にはGitHubのダウンロードページであるが，このオレンジでマークしたtar.gzファイルをダウンロードした方が早かったであろう。これは7zで解凍する。Windows64用である。

　図5〜7は，新たにダウンロードしたCloudCompare v. 2.13βの表情である。CloudCompareを選択して，v. 2.12.4と比べてみた。確実に改善されている。灰色だったものがactiveになっている。図５のように，開いたPluginsの下方にCloud layersが見える。これを選択した結果が，図６で，Cloud has no scaler fieldsのエラーが現れる。図７には，Editのメニューが見えている。このScaler fieldsを選択した際のリストが見えているが，このうち，Add constant SFなどが関係しているようである。このSFは，Scaler fieldsの略号と思われ，この定義をした上で，Cloud layersを見ることができるのだろう。

　このCloudCompareのCloud layersはこの程度のものなのか。とにかく，Danielの意見に戻らないといけない。

4. gcpクラウドを作成すべく既存のツールに期待したのだが

　次の，Edit > Cloudのうち，single point cloudしか，アクティブになっていない。もしクリップボードにデータがあれば，アクティブになる可能性がある。

• ‘Edit > Cloud > Create single point cloud‘

  • to create a cloud with a single point (set by the user)

• ‘Edit > Cloud > Paste from clipboard‘ (shortcut: CTRL+P)

  • to create a cloud from ASCII/text data stored in the clipboard

２章の末尾のDanielさんのアドバイスを次に再掲する。

Maybe you can use the new ‘qCloudLayers’ plugin?
Another way, maybe more simpler, is to use the ‘Edit > Scalar fields > Split cloud (integer values)‘ method. This will create as many (sub) clouds as classes present in your cloud.

　メニューを見ると，Edit > Scalar fields > Split cloud (integer values)，はグレイになっている。これも何らかの前処理があってのことである。さて。

以上，19:32，2024/01/07記。

　ちょっと，考えて見て，gcpのクラウドをエクセルやエディターで作成して，CloudCompareに取り込んで，メーンのdroneから得られたcloudの上のレイヤーとして表示するだけのことじゃないか，と考えついた。以前から頭にあったことが出てきたという感じである。

　で，その構造であるが，easting, northing, altitude,なのだけど，できれば，gcpの名称も欲しい。その観点で次に。

5. gcp cloudの書式は？

　どうしてこんな簡単なことをすぐにしなかったのか。CloudCompareでは幾つものメッシュやクラウドを取り込んで，DB Treeで，上下に移動してきたので，CloudCompareを使い始めてから，実はクラウドのレイヤーを操作してきたのである。

5.1 まずはEdit > Cloud > Create single point cloud

　結局，数点のクラウドをどのように作るかである。Edit > Cloud > Create single point cloud，って，何だろうか，ってわけで実際に使ってみた。x（easting: 地球座標系とは異なることに注意）, y（northing）, z値（altitude）を入力した（図９）。この値は，苦肉の策の，座標回転後のものである。

　この本題とは異なるが，2020年の調査結果をまとめたフォルダー毎，みつからない。別氏から送られてきたこれまでのdroneデータを全部削除する過程で，このフォルダーも殺したらしい。macのトラブルがあって，それに関するページを次のようにDec. 14に作成している。その翌日にTimeMachineでバックアップしたのだが，ここに殺したフォルダが残っていて，コピペした。幸い，座標回転までのエクセルデータがあった。このあたりから，別氏の再計算を待つ形で，キンドル本発行に注力したらしい。助かった。インフルエンザで多少熱があて鼻水だらだらながらの削除作業だったので，失策をしたらしい。怖いなあ−。

macのWi-Fiボタンが押せない　cannot push a button of wifi on mac

　図10は，座標値の簡略化の誘いの画面であるが，これに乗ってはならない。個々の座標値入力の際には必ず出てくる。drone測量の結果であるPLYファイルをCloudCompareに取り込む際にもNoを選択しているので，地球座標系の数値（GIS用にXYは例によって逆転）のままである。同じ座標系でないとgcpは宇和美崎オブジェクトに載らない。

　図10左ペーン最上部のDB Treeには，a_rot, f_rot，が見える。cloud☁︎（図中は白色だが）アイコンがついている。階層は🗂️アイコンがついたdroneファイルと同じだ。これでは，全部選んだとしても，CloudCompareのbinファイルとしては保存できない。

　ここで取り上げているCreate single point cloudコマンドで実行すると，DB Treeには自動で#001, #002，などとナンバーリングされる。この名称をタップすると，名称変更が可能なように入力できるようになって，a_rot，f_rot，とした訳である。メーンペーンでは白い小さな丸が現れて位置を確認できるのである。左ペーン中央のPropertiesのWindowでは，Name: a_rot, Visable: ☑️，Show Name：◻︎，が見える。Show Nameに☑️すると，メーンペーンに「a_rot」が現れる。

　これはいい！，と他の点も作成してゆくのである。5点まで登録して，「待てよ」となった。とにかく，一度，保存してみようと。図11左ペーンに見えるように，差し当たり作成した5点をまとめて保存することにした。File > save，だ。メーンペーンに，Save Ascii file，のパネルが現れた。しめしめ。パネルの意味はなんとなくわかる。

　図12　Save Ascii fileパネルに適宜入力したけど。coordinate precisionを8から12に（xy座標の期待できる精度から）。Scalarは変更せず6で（a_rotなどの名称ではないかと考えては英数字6文字で十分だと）。テキストとしては扱いやすいcsvで（それでseparatorをcommaに），白いアイコンも欲しいのでColorsもそのまま。

　さて，もともとのPLYがあったフォルダーには，5点が別々のcsvファイルとして保存されていた。名称は，gcp_000000.csv〜_000004.csv，である。まあ，いわばブラウザーでまとめて選んだだけで一つのファイルになるという想定が問題ではあった。

　なんと，ぼくはこのブラウザーで，gcp_000000.csvをgcp_a_rot.csvなどと変えちゃったりしたのである。これらのファイルを開くと自動でエクセルが開くが，名称やアイコンの色情報もなく，ただただ，x,y,zの値のみだ。

　CloudCompareで，図15のように，この5点をまとめて，読み込んだ。

　図16のように，パネル「Open Ascii file」が個々の点について現れるのだ。

　図17のように，5点がそれぞれ，フォルダーを持っていて，その下位に配置されている。一点ごとに。これは不合理だ。ポイント位置が見えない。そこで，Show nameに☑️すると，ポイントの名称が所定の場所に現れるのでおよそは位置がわかる。これではどうしようもない。PropertiesのBox dimensionsでは，x,y,zともゼロになっている。ここの「10」程度の数値が必要だろう。

　以上の経験から次の方針：　ポイント位置を保持するようなポイントの定義が必要だ。点をまとめて一つのフォルダーにまとめる必要がある。csvファイルに，ポイント名，座標情報，アイコン定義，が必要だ。

5.2 必要な属性を備えたgcpファイルを作成するには

　Create a new layer of just ground points (already classified)　での質問回答は，Danielが開発者故のものであった。フランスの質問者は，まさにぼくのこのテーマなのである。ベーシックな質問であった。

　ASCII export variables – what are they?　この情報に先に会っていたら，このページのここまでの試行錯誤は不要であったろう。図12のSave Ascii Fileパネルのパラメータに対するDanielの回答が次に。

Re: ASCII export variables – what are they?

Post by daniel » Sat Nov 05, 2022 3:03 pm

– coordinates precision: number of digits for coordinate values (X, Y, Z)
– scalar precision: number of digits for scalar field values (distance, intensity, etc.)
– separator: separator character between the above values (use a comma to create a CSV file for instance)
– order: order of the values (columns)
– header option: to add a description header (starts with ‘//’ and then you’ll find the name of each column)
– number of points: simply the total number of points on the first line (that’s expected if you want to create a valid ‘PTS’ file for some other tools)
– colors:
* use float values (between 0 and 1) instead of integer (between 0 and 255) for R, G and B components
* save the alpha (transparency) component as well

しかし，これではまだわからない。

　種々調べたけど，わからない。結局，自分でやるしかない。Danielさんに一度聞いてみたいものではあるが。xyz座標値は必須だが，マークbodyにxyz，それに，その色だ。そんなのはDanielさんの上のリストにはない。とにかく，エクセルでcsvファイルを作って，CloudCompareに読み込ませて，そこで指定してゆくのだろうなあ。

　図18のように，マックのエクセルでcsvファイルを作ってウィンドウズの所定のフォルダーに入れた。それをCloudCompareで開いたら，図19のように，念願の表Open Ascii Fileが出た。これが欲しかった。やってみないとわからない。図19の葡萄茶色の行はぼくが作成した列タイトルだ。適当に，Placemark, x-easting, y-northing, z-altitude, body-x, -y, -z, R, G, Bとしたのが，このテキストから適当に解釈して，それぞれ， この葡萄茶色の行のすぐ上の灰色の行のように，Label, cood.X, Y, Z, SF Scalar 3個，Red (0-255), Green (0-255), B (0-255), と配列している。

　自ら設定した列タイトルに合わない場合，図22のように，変更することができる。

　とにかく，図22の右下の，apply allのボタンを押した。図23の右上には正方形が表示されている。box-x,-y,-zとして全部10としたものである。図23で画像の拡大縮小に関わらずこの正方形のサイズは変わらない。アホだな。図24のように全7点が点アイコンもこの画面では変色しているが，指定の黄色だ。

　結局のところ，Scalar fieldの役割はわからないままだが，この3列を廃棄すればいいことがわかった。

6. drone cloudとgcpとの対応関係を見る

　上記のように, csvファイルの作成要素はわかったので，マックで編集して，図25のようにCloudCompareで読み込んだが，その結果が図26である。この画面では黄色というよりも明るい緑色になっている。Google Earthでもよく体験することだ。

　これ以降については，研究に流れ込むので，プライベートページに移動する。

以上，22:11，2024/01/08記。

1. 期待していること

　GrassGISをかつては使ってきたがご無沙汰している。Qgisは１回だけ使って精度面で問題を多く感じてアプリも削除した。Qgisは大衆化が著しく，悪貨が良貨を駆逐する，勢いである。日本の行政はGrassGISを飛び越して，専らQGIS派である。行政のコンテンツがQgisを意識したものになっている。

　さて，3Dメッシュまたは点群ファイルへの対応はどうだろうか。もちろん，GrassGISではこのアプリ公開以来，早くから米軍やNASAのLiDARデータの処理が実行されてきたのであるが，ぼくはこのコンテンツには触れなかった。

　いま，求めるのは，3Dモデルの3D表示能力はもちろんである。CloudCompareはその点では問題はない。しかし，CloudCompareでは，3Dモデルに，他のベクトルデータをオーバーレイすることができない。これはかなり不便なことである。それで，GISソフトに期待してのこのページ作成である。

　繰り返すと，3Dモデル表示が自在で，他のベクトルデータをオーバーレイできるか。この観点で，両アプリを調べたいと思う。

　

2. GrassGISではどうか

2.1 ケース１

　まずは次の回答を紹介する。

Import Point-Cloud to Grass-Gis
Data Processing – Discussion and Q&A

MichaelL

Feb ’20

Got it. I know Lidartools (which LAStools is dependent on) was not ported for QGIS 3 so that what I was getting at. 10M points isn’t all that big so if you have at least 16GB of RAM you shouldn’t have any problems, but considering you are having problems with a 500 point subset then we know that combined with the scaling issue that it’s probably not your plug-ins. Have you tried downloading and importing an XYZ version?

You might try PDAL.
https://grass.osgeo.org/grass78/manuals/v.in.pdal.html 23

　QGIS 3では対応できない。GLASSではXYZ形式であれば扱えるのではということである。RGB情報を破棄してしまうのは余りに不便でGLASSでも使い物にならないということだ。

2.2 r.in.pdal – Creates a raster map from LAS LiDAR points using univariate statistics

　このタイトルのリンクは，GLASSのマニュアルページである。

The r.in.pdal module loads PDAL library supported point clouds (with emphasis on LiDAR LAS files) into a new raster map using binning. The user may choose from a variety of statistical methods which will be used for binning when creating the new raster.

と，あって，精度が落ちるという点で，この種の処理をぼくは求めていないので，脚下だ。なお，ここでのテーマはLiDAR測量の生データの処理であり，いま，agisoftで，plyファイルなどが手に入るので，このケース１は求めるものではない。

2.3 ケース２

　それでは，まずは，plyファイルをGrassGISが取り込めるかについて調べてみる。

v.in.ply – Creates a vector map from a PLY file.

DESCRIPTION

v.in.ply imports a vector map in PLY vector format. A PLY file always holds a number of vertices which are imported as points. PLY vertices can have a number of properties in addition to their coordinates. These properties are stored in an attribute table. For larger PLY files with many vertices (> 1000) it is highly recommended to not use DBF as database driver, but SQLite (default in GRASS GIS 7), PostgreSQL or MySQL, because the DBF driver is rather slow and can consume a lot of memory. The database driver can be set with db.connect.

NOTES

v.in.ply is designed for large point clouds with the possibility to have only coordinates, and no attribute table (for speed reasons).

と，あって，興味深いことを次に。

　GISデータの基本コンテンツと考えていたDBFが，GrassGIS7以降，デフォルトとして，使われていないことである。SQLite (default in GRASS GIS 7)。これには驚いた。ぼくは時代遅れになっていた。
　上記Descriptionを見て大変期待できると思ったが，上記Notesを見ると，only coordinates, and no attribute table (for speed reasons)，とあって，ショックではあった。PLYファイルからRGB情報が削除されるのである。

　こう考えると，3D表示の可能性も見えないし，GISを使うメリットが感じられなくなったのである。

　v.ply.rectify – Imports PLY points, georeferences and exports them.　も重要だ。

DESCRIPTION

v.ply.rectify imports a PLY point cloud, georeferences and exports it. The first three vertex properties must be the x, y, z coordinates with property names “x”, “y”, “z”, in this order.

A text file with Ground Control Points (GCPs) must exist in the same folder where the point cloud is located, and the textfile must have the same name like the point cloud, but ending on .txt instead of .ply.

The text file with GCPs must have the following format with one GCP per line: x y z east north height status with x, y, z as source coordinates and east, north, height as target coordinates. The status indicates whether to use a GCP (status is not zero) or not (status is zero). Entries must be separated by whitespace or tabs. Decimal delimiters must be points.

The georecitifictation method used is a 3D orthogonal rectification where angles are preserved. 3D objects are shifted, scaled and rotated, but no shear is introduced. Please read the output of the module, in particular the root mean square (RMS) errors.

v.ply.rectify optionally exports the georeferenced point cloud not only with real coordinates, but also with shifted coordinates (-s flag) for display in meshlab or similar software that can not deal with real coordinates. The exported PLY point clouds will be in the same folder like the input PLY point cloud.

　次の観点は重要である：　x y z east north height status with x, y, z as source coordinates and east, north, height as target coordinates。通常のGISの座標系と異なり，x,y,z=easting, northing, heightになっている。別氏によればメタシェープでの処理も同様のようだ。

3. Qgisではどうか

　QGISで色付きの3Dメッシュを作成する は参考になったが，かなり，面倒な過程を経ている。わざわざ，他のカラー付きの画像を利用するのも理解できない。ここでは結局のところ，MeshLabが使われている。Qgisには3D表示能力がないようだ。この図１の右手のペーンには，レイヤーが配置されている。複数の3Dレイヤーを重ねて表示できるのである。不勉強でこの基礎的な機能についてさえ，知らなかった。

　Blenderでも可能である。Blenderで点群(.ply)を色付きでレンダリングする　Sep.23, 2023
出来上がった3Dオブジェクトを見ると，いかにも寂しい。

　で，SKラボ.netの，QGISで色付きの3Dメッシュを作成する，に出会った，と思った。しかし，元々RGB情報を持つplyを直接使わないのか，わからない。plyをインポートする際にRGB情報を廃棄してしまうからではあるが。結局，ぼくには無用のものであった。

4. MeshLabで

　結局，GLASSもQgisも使えないことがわかった。MeshLabに期待したい。MeshLabの操作法を学んで，このページに示したいと思う。

以上，19:15，2024/01/06記。

　Filters/Sampling/Clustered Vertex Subsampling, で，Cell Size1.0を例えば0.3などに変更してApplyしようとしたが，メーンメニューのFiltersをタップするだけで，応答無し，でフリーズしてしまう。重すぎるのだ。

　仕方無く，CloudCompareでFiles > Openを実行した。

5. CloudCompareでdown-samplingを

　試行錯誤してやっと見つけた。すごく簡単ではあった。Edit > subsample を実行すれば良い。開いたパネルで，min space btw points =，入力する数値であるが，2.7mにすると，出来上がったオブジェクトは7073点でスカスカだった。0.01mにすると全く変化なし。同じだ。で，0.100（10cm）にすると（図4は実行開始時, 図5は完成後），3,135,591点で，適当であった。元々２億点（198,925,288）からかなりスマートになった。

　図６には保存した新たなPLYファイルが見えていて，5GBから238MBへとよほど軽くなっている。

　では，MeshLabへ，という流れになるはずではあったが，偶然，Danielの凄いアドバイスを見つけてしまった。

Create a new layer of just ground points (already classified)　である。

以上，0:15，2024/01/07記。

6. CloudCompareでgcpも表示することに

　上記の，Create a new layer of just ground points (already classified) ,の内容を次にコピペする。

Create a new layer of just ground points (already classified)
Joined: Tue Aug 09, 2022 10:34 am
Create a new layer of just ground points (already classified)
Post by okay » Tue Aug 09, 2022 10:36 am

Hi -　I’m trying to make a separate layer of just the ground points in my point cloud. In other words, I want to be able to turn on and off a plot of just the ground points. I know about the CSF plug, but this isn’t necessary to use since the point cloud is classified. I can only figure out to how filter by value, but not by classification. Thanks.
Location: Grenoble, France
————————————————
Contact: Contact daniel
Re: Create a new layer of just ground points (already classified)
Post by daniel » Wed Aug 10, 2022 5:56 pm
Maybe you can use the new ‘qCloudLayers’ plugin?
Another way, maybe more simpler, is to use the ‘Edit > Scalar fields > Split cloud (integer values)’ method. This will create as many (sub) clouds as classes present in your cloud.
Daniel, CloudCompare admin
————————————————

　CloudCompareの開発者Danielさん，に，感謝。

　この情報を契機に，CloudCompareのPlugins Cloud layers，の存在を知った。図７のように，いま，使用中のCloudCompareの最上部には，CloudCompare v.2.12.4 (Kyiv) Stereo [64-bit]，が表示されている。メーンメニューPluginsに見えるCloud layersが欲しいのだが，灰色だ。

　このページのテーマが遷ってきているので，別にページを立てることにする。後に，そのページのリンクを示す。

以上，13:28，2024/01/07記。

CloudCompareでCloud layersを　handling cloud layers in CloudCompare