#!/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.1723632978996 56.7997390968573 -gcp 15171.0 0.0 25.3070839325154 56.8282053086749 -gcp 15171.0 17463.0 25.3709514407688 56.7528674337028 -gcp 0.0 17463.0 25.236230806153 56.7244012218852 LVVA_F7404_US1_GV2101_0002.jpg LVVA_F7404_US1_GV2101_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"312.0 13782.0\n516.0 12936.0\n624.0 12606.0\n996.0 11442.0\n1506.0 9294.0\n2292.0 6420.0\n2976.0 3648.0\n5136.0 2118.0\n6588.0 1308.0\n8850.0 510.0\n9300.0 1278.0\n10134.0 3126.0\n10854.0 4851.0\n11154.0 4791.0\n11793.0 5373.0\n12120.0 5559.0\n12342.0 5577.0\n12672.0 5442.0\n13035.0 5232.0\n13617.0 5004.0\n14394.0 4467.0\n14679.0 4257.0\n15165.0 3972.0\n15168.0 9408.75\n14145.0 10488.0\n13764.0 10827.0\n13503.0 11172.0\n13281.0 11547.0\n12990.0 11895.0\n12558.0 12282.0\n11112.0 13800.0\n10146.0 14124.0\n7596.0 15078.0\n5508.0 16020.0\n3252.0 16584.0\n1086.0 16650.0\n882.0 16338.0\n390.0 14958.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV2101_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV2101_0002.vrt LVVA_F7404_US1_GV2101_0002.cut.vrt

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