#!/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 24.322443008423 57.758128360556 -gcp 15060.0 0.0 24.459171295166 57.768658572718 -gcp 15060.0 8859.0 24.468612670898 57.728535781656 -gcp 0.0 8859.0 24.33274269104 57.716297718821 LVVA_F7404_US1_GV1203_0002.jpg LVVA_F7404_US1_GV1203_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"981.0 711.0\n3354.0 724.5\n3400.5 1008.0\n3699.0 1426.5\n4123.5 1875.0\n4698.0 2067.0\n5520.0 2145.0\n5811.0 2412.0\n6180.0 2817.0\n6603.0 3009.0\n7059.0 3216.0\n7359.0 3318.0\n7842.0 3159.0\n7938.0 3159.0\n8046.0 3378.0\n8268.0 3597.0\n8349.0 3735.0\n8430.0 3927.0\n8595.0 4083.0\n8742.0 4284.0\n8871.0 4497.0\n8991.0 4890.0\n9012.0 4983.0\n10173.0 4983.0\n10209.0 7821.0\n921.0 7737.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1203_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1203_0002.vrt LVVA_F7404_US1_GV1203_0002.cut.vrt

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