from osgeo import gdal, ogr, osr
import sys
import numpy as np

TIF_OLD = r"C:\Users\timon\Resilience BV\4020 SCane ESA DEMO - Documenten\General\4020 SCDEMO Team\4020 TechnicalData\WP3\smartcane_v2\smartcane\laravel_app\storage\app\angata\comparing_tiff\2025-11-25_old.tif"
TIF_OPT = r"C:\Users\timon\Resilience BV\4020 SCane ESA DEMO - Documenten\General\4020 SCDEMO Team\4020 TechnicalData\WP3\smartcane_v2\smartcane\laravel_app\storage\app\angata\comparing_tiff\2025-11-25_opt.tif"
GEOJSON = r"C:\Users\timon\Resilience BV\4020 SCane ESA DEMO - Documenten\General\4020 SCDEMO Team\4020 TechnicalData\WP3\smartcane_v2\smartcane\laravel_app\storage\app\angata\Data\pivot.geojson"


def geom_to_pixel_window(gt, geom):
    # bbox of geometry
    env = geom.GetEnvelope()  # returns (minX, maxX, minY, maxY)
    # Note: env is (minX, maxX, minY, maxY) but we'll map to pixel coordinates
    minx, maxx, miny, maxy = env[0], env[1], env[2], env[3]
    originX = gt[0]
    originY = gt[3]
    px_w = gt[1]
    px_h = gt[5]
    # compute pixel offsets
    xoff = int((minx - originX) / px_w)
    yoff = int((originY - maxy) / (-px_h))
    xend = int((maxx - originX) / px_w)
    yend = int((originY - miny) / (-px_h))
    xoff = max(xoff,0)
    yoff = max(yoff,0)
    xsize = max(1, xend - xoff + 1)
    ysize = max(1, yend - yoff + 1)
    return xoff, yoff, xsize, ysize


def check_tif_for_geom(tifpath, geom):
    ds = gdal.Open(tifpath)
    if ds is None:
        return None
    gt = ds.GetGeoTransform()
    band1 = ds.GetRasterBand(1)
    xoff, yoff, xsize, ysize = geom_to_pixel_window(gt, geom)
    try:
        arr = band1.ReadAsArray(xoff, yoff, xsize, ysize).astype(np.float32)
    except Exception as e:
        return {'error': str(e)}
    # check for any finite > 0 values (reflectance)
    finite = np.isfinite(arr)
    valid_pixels = np.sum(finite & (arr > 0))
    total_pixels = arr.size
    return {'valid_pixels': int(valid_pixels), 'total_pixels': int(total_pixels)}


if __name__ == '__main__':
    ds = ogr.Open(GEOJSON)
    lyr = ds.GetLayer()
    old_ds = gdal.Open(TIF_OLD)
    opt_ds = gdal.Open(TIF_OPT)
    if old_ds is None or opt_ds is None:
        print('Error opening one of the tiffs')
        sys.exit(1)

    gt_old = old_ds.GetGeoTransform()
    gt_opt = opt_ds.GetGeoTransform()

    results = []
    for feat in lyr:
        fid = feat.GetFID()
        geom = feat.GetGeometryRef()
        # compute windows
        old_res = check_tif_for_geom(TIF_OLD, geom)
        opt_res = check_tif_for_geom(TIF_OPT, geom)
        results.append({'fid': fid, 'old': old_res, 'opt': opt_res})

    # print summary
    found_old = sum(1 for r in results if r['old'] and r['old'].get('valid_pixels',0) > 0)
    found_opt = sum(1 for r in results if r['opt'] and r['opt'].get('valid_pixels',0) > 0)
    print(f'Total polygons checked: {len(results)}')
    print(f'Polygons with valid pixels in OLD tif: {found_old}')
    print(f'Polygons with valid pixels in OPT  tif: {found_opt}')

    # list differences
    missing_in_old = [r['fid'] for r in results if r['opt'] and r['opt'].get('valid_pixels',0) > 0 and (not r['old'] or r['old'].get('valid_pixels',0)==0)]
    missing_in_opt = [r['fid'] for r in results if r['old'] and r['old'].get('valid_pixels',0) > 0 and (not r['opt'] or r['opt'].get('valid_pixels',0)==0)]

    print('\nFields present in OPT but missing in OLD (FID list):')
    print(missing_in_old)
    print('\nFields present in OLD but missing in OPT (FID list):')
    print(missing_in_opt)

    # print per-feature counts for the first 10 features
    for r in results[:10]:
        print(f"FID {r['fid']}: OLD {r['old']}  OPT {r['opt']}")