#!/usr/bin/env bash

f="Bengaluru"
range='2023-04-01T00:00:00Z/2023-05-31T23:59:59Z'
srcbbox="61.966,5.051,96.595,33.946"
where="country = 'IN' and region = 'IN-KA' and \"names.primary\" like 'Bengaluru' and is_land = 1"
cloudcover="3"
poly="aoi-${f}.gpkg"

# Get location poly from overture
gdal vector pipeline \
   ! read -i /vsis3/overturemaps-us-west-2/release/2026-06-17.0/'theme=divisions/type=division_area' \
   ! filter --bbox ${srcbbox} --where "${where}" \
   ! clip --bbox ${srcbbox} --bbox-crs EPSG:4326 \
   ! write $poly --output-layer aoi --overwrite

# get bbox from poly extent
bbox="$(gdal vector info $poly --of json|jq -r '.layers[].geometryFields[].extent|@csv')"

# get sentinel meta data constrained by bbox and date range
gdal vector convert \
   -i "https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a&datetime=${range}&bbox=${bbox}&limit=500" \
   -o ${f}.geojson --overwrite

# SQL template for querying meta data, creating a list of all assets intersecting above poly
sql="select group_concat('[ ! read /vsicurl/' || \"assets.%s.href\" || ' ! reproject --bbox %s --bbox-crs EPSG:4326 --output-crs EPSG:4326 ]', ' ') from %s where \"eo:cloud_cover\" < %s"

for band in $(echo "red green blue") ; do
   gdal raster pipeline ! stack --resolution highest $(gdal vector sql -i ${f}.geojson --dialect sqlite --lco HEADER=NO --sql "$(printf "$sql" "${band}" "$bbox" $f $cloudcover)" --of CSV -o /vsistdout/|tr -d '"') ! write -o ${band}.tif --co COMPRESS=DEFLATE --overwrite
done

gdal raster pipeline \
   ! stack \
      [ ! calc -i "r=red.tif" --calc="avg(r) * .1 + .0001" --ot Byte --nodata 255 --flatten ] \
      [ ! calc -i "r=green.tif" --calc="avg(r) * .1 + .0001" --ot Byte --nodata 255 --flatten ] \
      [ ! calc -i "r=blue.tif" --calc="avg(r) * .1 + .0001" --ot Byte --nodata 255 --flatten ] \
   ! clip --like $poly --allow-bbox-outside-source \
   ! write -o ${f}.tif --co COMPRESS=DEFLATE --overwrite 
   gdal raster convert -i ${f}.tif -o ${f}.png --overwrite