#!/usr/bin/env python3

# Convert all .asc rasters in ./input from EPSG:4326 to EPSG:3031 COGs.
#
# GDAL Python bindings must be installed and compatible with the
# installed GDAL binaries.

from pathlib import Path
from osgeo import gdal

gdal.UseExceptions()


# ==================================================
# CONFIG
# ==================================================

BASE_DIR = Path.cwd()
INPUT_DIR = BASE_DIR / "input"
OUTPUT_DIR = BASE_DIR / "output"

SOURCE_CRS = "EPSG:4326"
TARGET_CRS = "EPSG:3031"
NODATA = -9999


# ==================================================
# MAIN
# ==================================================

def main():
    OUTPUT_DIR.mkdir(exist_ok=True)

    asc_files = sorted(INPUT_DIR.glob("*.asc"))

    if not asc_files:
        raise FileNotFoundError(f"No .asc files found in {INPUT_DIR}")

    for input_path in asc_files:
        output_path = OUTPUT_DIR / f"{input_path.stem}.tif"

        gdal.Warp(
            str(output_path),
            str(input_path),
            format="COG",
            srcSRS=SOURCE_CRS,
            dstSRS=TARGET_CRS,
            resampleAlg="bilinear",
            dstNodata=NODATA,
            creationOptions=[
                "COMPRESS=DEFLATE",
                "PREDICTOR=2",
                "BIGTIFF=IF_SAFER",
            ],
        )

        print(f"Created: {output_path.name}")

    print("Done.")


if __name__ == "__main__":
    main()