import os
from osgeo import gdal
import numpy as np
import subprocess


path_wbt = R"C:\Temp\whitebox_tools.exe"
outdir = "C:\Temp"
path_DEM = R"C:\Temp\dem.sdat"
max_iter = 20


def load_raster(path_raster):
    driver = gdal.GetDriverByName("SAGA")
    ds = gdal.Open(path_raster)
    band = ds.GetRasterBand(1)
    arr = band.ReadAsArray()
    return arr

# BreachDepressionsLeastCost command line
cmd_raw = " ".join([
    path_wbt,
    "--run=BreachDepressionsLeastCost",
    f"--dem={path_DEM}",
    "--min_dist",
    ])


# Run first iteration which will be used as a reference
ii = 0
outBDLC_ori = os.path.join(outdir, f"dem_BDLC_{ii}.sdat")
print(cmd_raw + f" --output={outBDLC_ori}")
subprocess.run(f"{cmd_raw} --output={outBDLC_ori}")
arr_ori = load_raster(outBDLC_ori)

# Rerun breaching until differences with the reference are found
while ii < max_iter:
    ii += 1
    outBDLC = os.path.join(outdir, f"dem_BDLC_{ii}.sdat")
    subprocess.run(f"{cmd_raw} --output={outBDLC}")
    arr_ii = load_raster(outBDLC)
    nb_diff = np.sum(arr_ori != arr_ii)
    print(f"{ii}/{max_iter} - Pixel differences: {nb_diff}")

    if nb_diff:
        print("Differences found!")
        break
    else:
        os.remove(outBDLC)
        os.remove(outBDLC[:-4] + "sgrd")