# Compare ArcGIS vs QGIS Flow Direction
# Run this in QGIS Python Console

from osgeo import gdal
import numpy as np

# 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.tif'
# arc_path = r'C:/Users/acer/Downloads/Documents/GIS_Work/DEM_Data/ArcGIS Hydrological Analysis/flow_dir_/vat.adf'
# qgis_path = r'C:/Users/acer/Downloads/Documents/GIS_Work/DEM_Data/QGIS_Output/flow_dir_6.tif'

print("=" * 60)
print("FLOW DIRECTION COMPARISON")
print("=" * 60)

print(f"ArcGIS Flow Direction: {arc_path.split("/")[-1]}")
print(f"QGIS Flow Direction: {qgis_path.split("/")[-1]}")

# Read ArcGIS Flow Direction
arc_ds = gdal.Open(arc_path)
if arc_ds:
    arc_band = arc_ds.GetRasterBand(1)
    arc_data = arc_band.ReadAsArray().astype(np.float64)
    arc_nodata = arc_band.GetNoDataValue()

    # Replace nodata with NaN
    if arc_nodata is not None:
        arc_data[arc_data == arc_nodata] = np.nan
    arc_data[arc_data < -1e30] = np.nan

    print('\n=== ArcGIS Flow Direction ===')
    print(f'Shape: {arc_data.shape}')
    print(f'NoData value: {arc_nodata}')
    valid_arc = arc_data[~np.isnan(arc_data)]
    print(f'Valid cells: {len(valid_arc):,}')
    print(f'Min: {np.min(valid_arc)}, Max: {np.max(valid_arc)}')

    # Count each D8 value
    print('\nValue distribution:')
    for val in [0, 1, 2, 4, 8, 16, 32, 64, 128]:
        count = np.sum(arc_data == val)
        pct = 100 * count / len(valid_arc) if len(valid_arc) > 0 else 0
        print(f'  {val:3d}: {count:>10,} cells ({pct:5.2f}%)')
else:
    print('ERROR: Could not open ArcGIS file')
    print(f'Path: {arc_path}')

# Read QGIS Flow Direction
qgis_ds = gdal.Open(qgis_path)
if qgis_ds:
    qgis_band = qgis_ds.GetRasterBand(1)
    qgis_data = qgis_band.ReadAsArray().astype(np.float64)
    qgis_nodata = qgis_band.GetNoDataValue()

    # Replace nodata with NaN
    if qgis_nodata is not None:
        qgis_data[qgis_data == qgis_nodata] = np.nan
    qgis_data[qgis_data < -1e30] = np.nan

    print('\n=== QGIS Flow Direction ===')
    print(f'Shape: {qgis_data.shape}')
    print(f'NoData value: {qgis_nodata}')
    valid_qgis = qgis_data[~np.isnan(qgis_data)]
    print(f'Valid cells: {len(valid_qgis):,}')
    print(f'Min: {np.min(valid_qgis)}, Max: {np.max(valid_qgis)}')

    # Count each D8 value
    print('\nValue distribution:')
    for val in [0, 1, 2, 4, 8, 16, 32, 64, 128]:
        count = np.sum(qgis_data == val)
        pct = 100 * count / len(valid_qgis) if len(valid_qgis) > 0 else 0
        print(f'  {val:3d}: {count:>10,} cells ({pct:5.2f}%)')
else:
    print('ERROR: Could not open QGIS file')
    print(f'Path: {qgis_path}')

# Compare if both loaded and same size
if arc_ds and qgis_ds:
    if arc_data.shape == qgis_data.shape:
        print('\n' + '=' * 60)
        print('COMPARISON')
        print('=' * 60)

        # Create masks for valid cells
        both_valid = ~np.isnan(arc_data) & ~np.isnan(qgis_data)

        # Compare values
        matching = np.sum((arc_data == qgis_data) & both_valid)
        total_comparable = np.sum(both_valid)

        print(f'\nCells valid in both: {total_comparable:,}')
        print(
            f'Matching direction:  {
                matching:,} ({
                100 * matching / total_comparable:.2f}%)')
        print(f'Different direction: {total_comparable -
                                      matching:,} ({100 *
                                                    (total_comparable -
                                                     matching) /
                                                    total_comparable:.2f}%)')

        # Analyze differences
        diff_mask = (arc_data != qgis_data) & both_valid
        diff_count = np.sum(diff_mask)

        if diff_count > 0:
            print('\nDifference Analysis:')
            # Sample some differences
            diff_rows, diff_cols = np.where(diff_mask)
            sample_size = min(10, len(diff_rows))
            print(f'\nSample differences (first {sample_size}):')
            print('  Row    Col   ArcGIS  QGIS')
            for i in range(sample_size):
                r, c = diff_rows[i], diff_cols[i]
                print(
                    f'  {r:5d}  {c:5d}  {int(arc_data[r, c]):5d}  {int(qgis_data[r, c]):5d}')
    else:
        print('\nERROR: Rasters have different shapes!')
        print(f'ArcGIS: {arc_data.shape}')
        print(f'QGIS: {qgis_data.shape}')

print('\n' + '=' * 60)
print('DONE')
print('=' * 60)
