89 lines
3.6 KiB
Python
89 lines
3.6 KiB
Python
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']}")
|