#!/bin/bash

GDAL_TILE_PROCESSES=16
GDAL_TILE_ZOOMS=8-14
GDAL_SAMPLING_WARP=cubic
GDAL_SAMPLING_TILE=antialias

# Create virtual dataset with coordinates
gdal_translate -of VRT -a_srs EPSG:4326 -outsize 300% 300% -gcp 0.0 0.0 25.688652992249 56.649503057559 -gcp 15060.0 0.0 25.830316543579 56.661651550098 -gcp 15060.0 9744.0 25.845551490784 56.613081318225 -gcp 0.0 9744.0 25.703887939454 56.600932825686 LVVA_F7404_US1_GV1675_0002.jpg LVVA_F7404_US1_GV1675_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"189.0 1311.0\n9480.0 1380.0\n9501.0 1491.0\n14602.5 1507.5\n14629.5 8493.0\n10116.0 8443.5\n10113.0 8350.5\n9504.0 8349.0\n6864.75 8334.75\n7386.0 8888.25\n5748.0 8886.75\n5663.25 8329.5\n183.0 8349.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1675_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1675_0002.vrt LVVA_F7404_US1_GV1675_0002.cut.vrt

# Generate tiles
gdal2tiles.py -r $GDAL_SAMPLING_TILE --xyz -z $GDAL_TILE_ZOOMS -x --processes=$GDAL_TILE_PROCESSES LVVA_F7404_US1_GV1675_0002.cut.vrt LVVA_F7404_US1_GV1675_0002.xyz