diff --git a/python_app/02_harvest_imminent_weekly.py b/python_app/31_harvest_imminent_weekly.py similarity index 100% rename from python_app/02_harvest_imminent_weekly.py rename to python_app/31_harvest_imminent_weekly.py diff --git a/tools/compare_tifs.py b/tools/compare_tifs.py deleted file mode 100644 index 5eec4e1..0000000 --- a/tools/compare_tifs.py +++ /dev/null @@ -1,88 +0,0 @@ -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']}")