#!/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.663721621669 57.322476216233 -gcp 13077.0 0.0 25.755043029785 57.274356509677 -gcp 13077.0 9846.0 25.692043304443 57.239678255543 -gcp 0.0 9846.0 25.600110225169 57.285949316601 LVVA_F7404_US1_GV2008_0002.jpg LVVA_F7404_US1_GV2008_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"429.0 3655.5\n796.5 2557.5\n1165.5 1641.0\n1327.5 1581.0\n1338.0 1257.0\n1398.0 1194.0\n2880.0 1021.5\n3156.0 955.5\n3160.5 600.0\n3433.5 490.5\n3727.5 468.0\n3901.5 175.5\n4119.0 180.0\n4321.5 238.5\n4554.0 597.0\n4816.5 568.5\n5041.5 604.5\n5394.0 732.0\n6109.5 666.0\n6852.0 615.0\n7288.5 562.5\n7785.0 459.0\n8604.0 360.0\n9555.0 1201.5\n9895.5 777.0\n10000.5 744.0\n10191.0 433.5\n10755.0 501.0\n10998.0 603.0\n10978.5 831.0\n12090.0 1342.5\n12012.0 2433.0\n12036.0 3379.5\n12111.0 4455.0\n12225.0 4938.0\n12525.0 5718.0\n12864.0 6420.0\n12873.0 6567.0\n12684.0 7180.5\n12490.5 7803.0\n11542.5 8412.0\n9996.0 8937.0\n7929.0 8925.0\n6792.0 8619.0\n6141.0 8169.0\n5727.0 7686.0\n5406.0 7203.0\n4098.0 7128.0\n3246.0 7287.0\n2490.0 7326.0\n2014.5 7297.5\n1650.0 7197.0\n520.5 6559.5\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV2008_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV2008_0002.vrt LVVA_F7404_US1_GV2008_0002.cut.vrt

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