Three-Dimensional Back Analysis コード修正

文献：Geologic Trends in Shear Strength Properties Inferred Through Three‐Dimensional Back Analysis of Landslide Inventories - Bunn - 2020 - Journal of Geophysical Research: Earth Surface - Wiley Online Library

実装：GitHub - benalesh/Landslide-Forensics の修正箇所

BackAnalysis_3D.m

% Read attributes from F
％ポリゴンに属性を与えていない場合の代替え処理
%asp = F(idx).aspect;% Mean aspect of all surficial pixels
%fail_type = F(idx).FAIL_TYPE;% Fail type: 'Catastrophic' or 'Progressive'
asp = 0; %すべり方向は後に繰り返し処理で決定する
fail_type = 'Progressive'

% Create binary mask from input polygons
%[x,y] = pixcenters(geotiffinfo);%現在は廃止された関数
%[X,Y] = meshgrid(x,y);
[X,Y] = worldGrid(R);

% Loop through all features in F (extent shapefile)
for idx = 1:length(F)
    try %エラー処理追加_
        繰り返し処理
  catch ME
  warning('Slide %d: Error occurred. Skipping...\n\t%s', idx, ME.message);
        % 必要に応じて、失敗スライド用の初期値を入れておく (任意)
        F(idx).cf = double(0);
        F(idx).phif = double(0);
        F(idx).rot = double(0);
     continue;
    end
end

% F 内のすべてのフィールド名を取得
fn = fieldnames(F);
% F の各属性に対してループ
for iF = 1:numel(F)
    for iField = 1:numel(fn)
        fieldName = fn{iField};
        val = F(iF).(fieldName);
        % 数値型かつ空配列なら NaN で埋める
        if isnumeric(val) && isempty(val)
            F(iF).(fieldName) = NaN;
        end
    end
end

%ヒストグラムの作成
phif_arr=[F.phif]
idx_valid  = (phif_arr > phi_thresh); % 論理索引（true/false配列）
phif_filt  = phif_arr(idx_valid); % 有効な値のみ抽出
figure;
h = histogram(phif_filt, 10, 'Normalization','probability');
xlabel('phif (deg)');
ylabel('Probability Density');
title('Distribution of phi');
% h.Values : 各ビンの高さ (確率)
% h.BinEdges : ビン境界(N+1個)
% ビンの中心を各隣り合う境界の平均で計算
center_phi = (h.BinEdges(1:end-1) + h.BinEdges(2:end)) / 2;
prob     = h.Values;
prob_phi = center_phi;
prob_coh  = repmat(c_thresh, 1, 10);
save(fullfile(path, 'shear_strength.mat'), 'prob', 'prob_coh', 'prob_phi');

SimpJanbu3D.m

    maxC   = 1000;     % c の上限 [kN/m^2]
    maxPhi = 90;       % φ の上限 [度]
    while FRy<FDy
        if st == 1
            phi = phi+phi_inc;
            if phi >= maxPhi
                break;  % 上限到達で打ち切り
            end
        else
            c = c+c_inc;
            if c >= maxC
                break;  % 上限到達で打ち切り
            end
        end