#!/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.176372528076 57.542733466419 -gcp 14265.0 0.0 25.31533241272 57.559588548886 -gcp 14265.0 14067.0 25.343570709229 57.495350142975 -gcp 0.0 14067.0 25.210018157959 57.476065775889 LVVA_F7404_US1_GV626_0002.jpg LVVA_F7404_US1_GV626_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"2508.0 8235.0\n699.0 8241.0\n693.0 4128.0\n447.0 3468.0\n441.0 2992.5\n705.0 2838.0\n711.0 2091.0\n2427.0 2058.0\n2436.0 264.0\n11442.0 417.0\n11442.0 873.0\n11650.5 879.0\n11947.5 565.5\n12903.0 589.5\n13141.5 868.5\n13743.0 882.0\n13698.0 6987.0\n12109.5 7000.5\n11757.0 7503.0\n11749.5 7767.0\n11458.5 7929.0\n11463.0 8724.0\n11752.5 8832.0\n11733.0 9955.5\n11457.0 10020.0\n11472.0 13164.0\n2538.0 13101.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV626_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV626_0002.vrt LVVA_F7404_US1_GV626_0002.cut.vrt

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