地形の曲率

西田ほか「数値地形モデルに基づく地震時山腹崩壊斜面の地形解析」1997

いくつかの国総研資料にて、この文献による「平均曲率」が利用されています。文献でいう４方向ラプラシアンでは斜度により数値が変化するため、変化しない平均曲率が良いとのことでした。

式自体は汎用ソフトに実装されている式と似ています。が、５mDEMの利用で計算メッシュの距離を10m～50m とすると、カーネルが 3x3 では収まりません。これ、GISソフトで対応しているものがないのでは？

まず、ArcGIS。この文献で否定されている式のみの実装なのでダメですね。
SAGA には複数の曲率計算方法が実装されていますが、該当するものはないようです。
https://sourceforge.net/p/saga-gis/code/ci/master/tree/saga-gis/src/tools/terrain_analysis/ta_morphometry/Morphometry.cpp
GISではないものの、Surfer にもいくつか曲率の計算が実装されています。ウェーブレットで９x９のユーザー定義フィルターを使えましたので、もしかするとできるのかもしれません。が、手元になかったので確認できませんでした。
http://surferhelp.goldensoftware.com/gridops/profile_curvature.htm

数字の並びにフィルタをかけていくのは画像処理と同じ作業です。OpenCV であれば、どのような間隔・カーネルでも対応可能です（しかも計算が速い）。いつもの Python だと、手順は以下の通りです。

1. GDAL で GeoTIFF 読み込み
2. NumPy でカーネル作成
3. OpenCV で ２D フィルタをかけて ∂h/∂x ほかを計算
4. H を計算
5. GDAL で GeoTIFF 書き出し

簡単、と思って手を動かし始めたのですが、詰まりました。
いえ、すべて正しく組めて速く走るのですが、できた GeoTIFF を ArcGIS や QGIS に読み込ませると最大・最小の表示が狂います（値自体は正しく読み込まれています）。
一晩考えたものの、原因を特定できませんでした。

今朝、同じ問題で悩まれていた方を発見。これは知りませんでした。
https://gis.stackexchange.com/questions/29064/how-to-get-gdal-to-create-statistics-for-gtiff-in-python

6. 統計量を XNL で出力

最後にこの手順を加えると、あっさり解決。
今まで Geopandas では対応できないラスターの読み書きを避けていましたので、付けが回ってきたのでしょう。

やはり、手を動かすと、何かしらの発見があるものです。

**********************
20220731追記

ArcGIS Pro の エクステンションに追加されたようです。
https://pro.arcgis.com/ja/pro-app/2.8/tool-reference/spatial-analyst/surface-parameters.htm

平均曲率は Minár(2020)とされていますが、その文献の参照先は Gauss（1828）になっていました。

**********************
20220802追記
ラスターでもフィルタの距離を指定できました。
計算結果は±を逆転しても一致しません。

「ローカルサーフェスタイプ」を指定して、サーフェスを作成してから、正規傾斜ライン？＝コンター直交方向と接線方向を求めXY軸に定めているようです。谷底ではどうなるのでしょうか。