5mメッシュ（標高）の精度 その2

ある図郭をSurferで図化しました（griddingは約5m）。

5mメッシュ（LP）

 5mメッシュ（写真測量）

 見た目はさほど変わりません。が、差分（写真標高-LP標高）をとると、以下の通り。

差分（写真-LP）

やはり、大きく差が出ています。特に山中の谷部で20m以上の差が出ています。ま、木があって難しいんでしょうね。測量屋さんによると、写真で分からないところは、ある程度現地で確認して補正すると言われていましたが、広域の測量時に一々、道のない谷に入って確認するとも思えません。尾根でも10m以上の差の出ている箇所は結構あるので、単純に植生密度の影響かもしれません。

gridding 前の生データ（といっても、これも gridding 済みですが）を比較しても、写真測量が平均値、中央値で+10m、最大+28mとなっています。基本的に、データの質として、それだけの差があるということです。

当然、技術上も経験上も LP の方が実測に近いと思われますので、写真測量の精度が粗いと考えるべきなのでしょうね。

提供されているデータは、あくまで標高点のみ精度が保証されているということでしょう。写真測量の他の図郭においても、20m程度の差は含まれる可能性があると考えておいた方が良いかも知れません。 

Wavelet その2

メキシカンハット関数を使った、ウェーブレット解析（ただの演算ですが）の続きです。

昨日UPした通り、(1-x^2・・・から始めることにしました。
2次元関数はこのような感じです。

マニュアルではスケール S を1で固定し、DEM のグリッドサイズで波長が変わるといった内容です。λ=4グリッドといったところですね。実務的ですね。
x-aが影響するところまでを考えると、最低でも中心から左右に4グリッドずつ必要です。従って、フィルタは9*9のマトリックスになります。

早速、SAGA で User Defined Filter があるか調べてみました。結果、装備されていました。しかし、3*3までなので適用できません。仕方ないですね。
Surfer では9*9でも可能ですから、こちらでやってみましょう。ただし、Filter による演算結果は、最後にマトリックスの合計値で割る設定になっています。そのため、マニュアルの積分に合わせるには、後で Grid Math で合計値を掛け戻してやる必要があります。また、Filter の edge の処理設定もしないといけません。

演算結果を確認してみますと、あれ？といった感じの仕上がりです。全くマニュの図とは異なっています。係数のオーダーも大きく違います。チェックしましたが計算過程は問題ありません。
ということは、(1-x^2・・・ではなくて、（2-x^2が正しいということでしょうか？

(2-x^2・・・の2次元関数は以下のようになります。振幅の最大値は大きくなりますが、同じような形状です。

Surfer で演算し直し、できた GRD ファイルを SAGA で読み込んだ結果がこちら。

オーダーもあっています。（2-X^2・・・で問題ないようでした。
最終的に、地すべり地形を判読しやすいと実感したら、式の導出まで追いかけてみましょうか。このまま使うのは気持ち悪いですからね。


とりあえず、地上開度、ウェーブレット解析、両方ともクリアーです。
地上開度もよく読んでみますと、1グリッドで計算した場合も書かれていますので、最大値を取るアルゴを組まなくてもよくなっています。非常に実務的と言いますか、汎用ソフトの利用を前提に、適用しやすいよう単純化して書かれています。技術的に良いか悪いかは分かりませんが、裾野を広げることはできるでしょう。


あとは地上開度図とウェーブレット解析図の合成のみです。
こちらはまた後日。

Wavelet

地上開度がクリアできたので、俄然やる気が出てきました。

土研「地すべり地における航空レーザー測量解析マニュアル（案）」に従えば、あと視覚化で残るのはウェーブレット解析です。

ウェーブレット関数とは何ぞや？の状態ですので、ネットで調べてもよく分かりませんでした。これだけで１冊の本になるほど多肢にわたっています。

しかし、マニュのメキシカンハット関数やウェーブレット係数の算出式は非常に簡素なので、これに限ってはできるような気がしました。

で、とりあえず、マニュの式をEXCELでグラフ化。

あれ？最大値が11ページのグラフと違います。おかしい。

ネットで調べてみると、ψ(x, y) = （2-x^2・・・ の箇所が（1-x^2・・・の形になっている式が多く見られます。y = 0のときに（1-x^2と同じ形にならないといけないはずですから誤植でしょうか？

とりあえず2を1に修正してみると、

11ページと同じ図ができました。私の知識ではどちらが正解かわかりませんが、係数を出すときの重みの大きさが変わるだけなので、とりあえずこちらで行きましょう。今度、地形解析チームに出典の文献を持っていないか聞いてみましょう。

あとは積分の式をどう反映させるかです。これはソフト側でユーザー定義フィルターが使用できるかどうかも重要でしょう。SAGA で可能なのでしょうか？Surfer では可能ですし、ImageJでもMexican Hat の計算はできますから、最悪はそちらでの検討になるでしょう。

続きは後日。

SAGA

SAGAの続きです。

LPデータを読み込みました。
ナンバリングしてあれば、SurferでXYZのみに加工したデータをCSVとして保存しておけば読み込めます(EXCEL2010では行数に限度があります)。また、Surferのグリッドデータ読めます。

Terrain Analysis Module がたくさん装備されています。
http://saga-gis.org/geostat2011/05_saga_geomorphometry.pdf
が、まだその計算内容のまとまった説明を見つけることができません。先日のようにWikiで一つ一つ探さないといけないのでしょうか？また、これらの計算は何の目的で考案されたのでしょうか？何がわかるのでしょうか？これをツールとして使う分野は？地形学とも違うような気がします。知りたいですね。

とりあえず、正常に動きましたので、作業自体は問題ないでしょう。


これから計算内容を確認していきましょう。