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