"""Synthetic raster exercise. Requires GDAL Python bindings and NumPy.
No drone imagery, field accuracy claim, or external network required.
Run: python gdal-raster-check.py
Temporary GeoTIFFs are removed when the script exits.
"""
import json
import tempfile
from pathlib import Path
import numpy as np
from osgeo import gdal, osr

gdal.UseExceptions()
checks = []


def check(name, condition):
    checks.append({"name": name, "passed": bool(condition)})
    if not condition:
        raise AssertionError(name)


with tempfile.TemporaryDirectory(prefix="drone-raster-") as scratch:
    source = str(Path(scratch) / "synthetic-utm.tif")
    target = str(Path(scratch) / "synthetic-wgs84.tif")
    crs = osr.SpatialReference()
    crs.ImportFromEPSG(32647)
    grid = np.full((10, 10), 10, dtype=np.float32)
    grid[0, 0] = -9999
    transform = (500000, 2, 0, 1500000, 0, -2)
    ds = gdal.GetDriverByName("GTiff").Create(source, 10, 10, 1, gdal.GDT_Float32)
    ds.SetSpatialRef(crs)
    ds.SetGeoTransform(transform)
    ds.GetRasterBand(1).SetNoDataValue(-9999)
    ds.GetRasterBand(1).WriteArray(grid)
    ds = None  # Close before reopening; works on GDAL 3.8 as well.
    ds = gdal.Open(source)
    check("dimensions_10_by_10", ds.RasterXSize == 10 and ds.RasterYSize == 10)
    check("source_epsg_32647", ds.GetSpatialRef().GetAuthorityCode(None) == "32647")
    check("geotransform_preserved", ds.GetGeoTransform() == transform)
    check("pixels_preserved", np.array_equal(ds.ReadAsArray(), grid))
    check("nodata_preserved", ds.GetRasterBand(1).GetNoDataValue() == -9999)
    valid = ds.GetRasterBand(1).GetMaskBand().ReadAsArray() != 0
    check("valid_cells_99", np.count_nonzero(valid) == 99)
    area = int(np.count_nonzero(valid)) * abs(transform[1] * transform[5])
    check("projected_valid_area_396_m2", area == 396)
    geographic = osr.SpatialReference()
    geographic.ImportFromEPSG(4326)
    for srs in [crs, geographic]:
        srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    lon, lat, _ = osr.CoordinateTransformation(crs, geographic).TransformPoint(500000, 0)
    check("utm_zone47_equator_reference", abs(lon - 99) < 1e-8 and abs(lat) < 1e-8)
    warped = gdal.Warp(target, ds, dstSRS="EPSG:4326", resampleAlg="near",
                       srcNodata=-9999, dstNodata=-9999, multithread=False)
    warped = None
    out = gdal.Open(target)
    check("target_epsg_4326", out.GetSpatialRef().GetAuthorityCode(None) == "4326")
    values = out.ReadAsArray()
    usable = values[values != -9999]
    check("nearest_preserves_constant_valid_values", usable.size > 0 and np.all(usable == 10))
    check("target_nodata_metadata", out.GetRasterBand(1).GetNoDataValue() == -9999)
    result = {"gdal_version": gdal.VersionInfo("RELEASE_NAME"), "fixture": "synthetic 10x10 constant raster with one NoData cell",
              "checks": checks, "source_valid_area_m2": area,
              "target_size": [out.RasterXSize, out.RasterYSize],
              "target_geotransform": out.GetGeoTransform(),
              "scope": "GeoTIFF and CRS software exercise only; not photogrammetry or physical validation"}
    out = None
    ds = None
print(json.dumps(result, indent=2))