import numpy as np
from osgeo import gdal

# File paths
arc_path = r'c:/Users/acer/Downloads/Documents/GIS_Work/DEM_Data/clip/clip_flow_dir/vat.adf'
qgis_path = r'C:/Users/acer/Downloads/Documents/GIS_Work/DEM_Data/clip/clipped_flow_dir_8.tif'
dem_path = r'C:/Users/acer/Downloads/Documents/GIS_Work/DEM_Data/Clipped_DEM/Clipped_DEM.0.tif'

print("=" * 60)
print("PIXEL-BY-PIXEL FLOW DEBUGGER")
print("=" * 60)
print(f"ArcGIS: {arc_path.split('/')[-1]}")
print(f"QGIS:   {qgis_path.split('/')[-1]}")
print(f"DEM:    {dem_path.split('/')[-1]}")
print("-" * 60)

# Open datasets
arc_ds = gdal.Open(arc_path)
qgis_ds = gdal.Open(qgis_path)

if not arc_ds or not qgis_ds:
    print("Error: Could not open one or more flow direction files.")
else:
    arc_band = arc_ds.GetRasterBand(1)
    qgis_band = qgis_ds.GetRasterBand(1)

    # Read arrays
    arc_arr = arc_band.ReadAsArray()
    qgis_arr = qgis_band.ReadAsArray()

    rows, cols = arc_arr.shape

    # Find mismatch locations
    # Ignore NoData (255 for ArcGIS usually, checking simplified logic)
    # Also ignore 0 if that's considered background
    diff_mask = (arc_arr != qgis_arr) & (arc_arr != 255) & (qgis_arr > 0)

    mismatch_indices = np.where(diff_mask)
    count = len(mismatch_indices[0])

    mismatch_indices = np.where(diff_mask)
    count = len(mismatch_indices[0])

    # Force full array printing (Global setting)
    np.set_printoptions(
        threshold=np.inf,
        linewidth=1000,
        edgeitems=1000,
        suppress=True)

    if count == 0:
        print("No mismatches found!")
    else:
        print(f"Total Mismatches: {count}")
        print("-" * 60)

        # Focus on the first mismatch
        center_r, center_c = mismatch_indices[0][0], mismatch_indices[1][0]

        # Define window size (Radius around center)
        # w_size = 2  -> 5x5 grid (2 left + center + 2 right)
        # w_size = 5  -> 11x11 grid
        w_size = 5
        r_start = max(0, center_r - w_size)
        r_end = min(rows, center_r + w_size + 1)
        c_start = max(0, center_c - w_size)
        c_end = min(cols, center_c + w_size + 1)

        print(f"Inspecting Window around ({center_r}, {center_c})")

        print("\n--- ArcGIS Flow Direction ---")
        print(arc_arr[r_start:r_end, c_start:c_end])

        print("\n--- QGIS Flow Direction ---")
        print(qgis_arr[r_start:r_end, c_start:c_end])

        # DEM Analysis
        if dem_path:
            dem_ds = gdal.Open(dem_path)
            if dem_ds:
                dem_arr = dem_ds.GetRasterBand(1).ReadAsArray()
                print("\n--- Elevation Values ---")
                subset = dem_arr[r_start:r_end, c_start:c_end]

                # Format for readability
                # np.set_printoptions(precision=3, suppress=True, linewidth=120)
                print(subset)

                # Analyze the center cell detailed math
                print(
                    f"\n--- Detailed Slope Analysis for ({center_r}, {center_c}) ---")
                curr_elev = dem_arr[center_r, center_c]
                print(f"Center Elev: {curr_elev:.4f}")

                drs = [0, 1, 1, 1, 0, -1, -1, -1]
                dcs = [1, 1, 0, -1, -1, -1, 0, 1]
                codes = [1, 2, 4, 8, 16, 32, 64, 128]
                names = ['E', 'SE', 'S', 'SW', 'W', 'NW', 'N', 'NE']

                geotransform = dem_ds.GetGeoTransform()
                cellsize_x = abs(geotransform[1])
                cellsize_y = abs(geotransform[5])
                # Simple assumption for metric projection usually
                # If projected, cellsize is in meters. calculate diagonal.
                diag = np.sqrt(cellsize_x**2 + cellsize_y**2)
                dists = [
                    cellsize_x,
                    diag,
                    cellsize_y,
                    diag,
                    cellsize_x,
                    diag,
                    cellsize_y,
                    diag]

                print(
                    f"{
                        'Dir':<4} {
                        'Code':<4} {
                        'Neighbor':<10} {
                        'Drop':<10} {
                        'Dist':<8} {
                            'Slope':<10}")
                print("-" * 60)

                max_s = -float('inf')
                best_code = 0

                print(
                    f"{
                        'Dir':<4} {
                        'Code':<4} {
                        'Neighbor':<10} {
                        'Drop':<10} {
                        'Dist':<8} {
                            'Slope':<10}")
                print("-" * 60)

                # We will collect all max_slope directions to see the tie set
                max_s = -float('inf')
                candidates = []

                for i in range(8):
                    nr, nc = center_r + drs[i], center_c + dcs[i]
                    if 0 <= nr < rows and 0 <= nc < cols:
                        n_elev = dem_arr[nr, nc]
                        drop = curr_elev - n_elev
                        slope = drop / dists[i]

                        is_max = False
                        if slope > max_s + 0.000001:  # Clear new max
                            max_s = slope
                            candidates = [codes[i]]
                            is_max = True
                        elif abs(slope - max_s) < 0.000001:  # Tie
                            candidates.append(codes[i])
                            is_max = True

                        marker = "*" if is_max else ""
                        print(
                            f"{
                                names[i]:<4} {
                                codes[i]:<4} {
                                n_elev:<10.4f} {
                                drop:<10.4f} {
                                dists[i]:<8.2f} {
                                slope:<10.6f} {marker}")

                print("-" * 60)
                print(f"Tie Candidates: {candidates}")
                print(f"ArcGIS Value:   {arc_arr[center_r, center_c]}")
                print(f"QGIS Value:     {qgis_arr[center_r, center_c]}")
                print(
                    f"Position: ({center_r}, {center_c}) Sum={
                        center_r +
                        center_c} Parity={
                        (
                            center_r +
                            center_c) %
                        2}")

                # Hypothesis Check
                # Does ArcGIS Pick Max or Min code based on parity?
                if arc_arr[center_r, center_c] in candidates:
                    chosen = arc_arr[center_r, center_c]
                    is_max_code = (chosen == max(candidates))
                    is_min_code = (chosen == min(candidates))
                    print(f"ArcGIS chose: {chosen}")
                    if is_max_code and is_min_code:
                        print("Method: Single candidate (No tie)")
                    elif is_max_code:
                        print("Method: MAX code (128 > ... > 1)")
                    elif is_min_code:
                        print("Method: MIN code (1 > ... > 128)")
                    else:
                        print("Method: Middle code (Complex priority?)")
                else:
                    print(
                        ">> WARNING: ArcGIS chosen value NOT in calculated candidates! Precision issue or different algorithm (Aspect?)")

                if best_code != arc_arr[center_r, center_c]:
                    print(">> MISMATCH: Logic differs from ArcGIS")
                else:
                    print(">> MATCH: Logic matches ArcGIS")

            else:
                print(f"Could not open DEM: {dem_path}")

print("\nDONE")
