3次スプライン関数で崩壊面推定 その2

Assessing the landslide failure surface depth and volume: A new spline interpolation method - ScienceDirect

GitHubにコードあり、と書かれていましたが、本日時点で ReadMe と License しか公開されていません。スプラインを使うだけの簡単な手法なので、書いてみました。

import numpy as np
from scipy.interpolate import CubicSpline
from tqdm import tqdm

def spline_failure_surface_pointwise(
    dem: np.ndarray,
    landslide_mask: np.ndarray,
    fix_points=None,
    max_iter=1000,
    tol=0.01,
    slope_angle_thresh=1
):

    failure_surface = dem.copy()

    fixed_mask = np.zeros_like(dem, dtype=bool)
    fixed_depths = np.zeros_like(dem)
    if fix_points:
        for (i, j), val in fix_points.items():
            fixed_mask[i, j] = True
            fixed_depths[i, j] = val

    nrows, ncols = dem.shape

    print('dem.shape:', dem.shape)
    print('landslide_mask.shape:', landslide_mask.shape)

    for it in tqdm(range(max_iter)):
        prev_surface = failure_surface.copy()

        for i in range(nrows):
            for j in range(ncols):
                if not landslide_mask[i,j]:
                    continue

                if fixed_mask[i,j]:
                    failure_surface[i,j] = fixed_depths[i,j]
                    continue

                j_idx_list = [jj for jj in range(ncols)
                             if landslide_mask[i,jj] and jj != j]
                if len(j_idx_list) >= 2:
                    row_vals = [failure_surface[i,jj] if not fixed_mask[i,jj] else fixed_depths[i,jj]
                                for jj in j_idx_list]
                    row_spline = CubicSpline(j_idx_list, row_vals, bc_type='not-a-knot')
                    zintp_i = row_spline(j)
                else:
                    zintp_i = failure_surface[i,j]

                i_idx_list = [ii for ii in range(nrows)
                             if landslide_mask[ii,j] and ii != i]
                if len(i_idx_list) >= 2:
                    col_vals = [failure_surface[ii,j] if not fixed_mask[ii,j] else fixed_depths[ii,j]
                                for ii in i_idx_list]
                    col_spline = CubicSpline(i_idx_list, col_vals, bc_type='not-a-knot')
                    zintp_j = col_spline(i)
                else:
                    zintp_j = failure_surface[i,j]

                # 両方向平均
                zavg = (zintp_i + zintp_j) / 2

                # 文献式:
                if zavg < failure_surface[i, j]:
                    failure_surface[i, j] = zavg

        # 収束判定（標高変化量）
        change = np.abs(failure_surface - prev_surface).max()
        if change < tol:
            print(f"Converged at iteration {it}")
            break

    return failure_surface

制御点をある関数にあてはめて残りを推定するというのはGEORAMA等と同じ考え方。GEORAMAはGUI（CAD）を利用しているので制御点を与えやすいのがメリット。

適当な地形を作って動かしてみたのが以下。

結果はイマイチ。イタレーションを3000回にしたのですが、まだ収束していません。これでも数時間かかったのですが、並列化しても半分ぐらいでしょう。DEMの解像度が高い、あるいは崩壊範囲が広いと時間のかかる手法だということがわかりました。こうなると、スプラインよりもGEORAMAの方が優秀、SLBLの方が簡素。関数を変更したら良いのでしょうね。https://www.blogger.com/blog/post/edit/269984068930159765/2223718948702813275

ま、ボーリングによる崩壊面深度を反映する場合の簡易手法ということで、覚えておきましょう。