#!/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.930124282837 57.66487101148 -gcp 9972.0 0.0 25.096549987793 57.660922682565 -gcp 9972.0 7296.0 25.098266601562 57.608681687188 -gcp 0.0 7296.0 24.92712020874 57.609923109825 LVVA_F7404_US1_GV269_0002.jpg LVVA_F7404_US1_GV269_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"556.5 2580.0\n273.0 2466.0\n241.5 1965.0\n358.5 1501.5\n399.0 1342.5\n525.0 1278.0\n510.0 363.0\n3295.5 358.5\n4723.5 355.5\n9534.0 384.0\n9490.5 4826.25\n9540.0 4848.75\n9558.0 4890.75\n9588.75 4914.75\n9489.0 4941.0\n9487.5 5037.0\n8472.0 5028.0\n8490.0 6626.25\n6637.5 6607.5\n6641.25 5376.0\n6228.0 5037.0\n576.0 5037.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV269_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV269_0002.vrt LVVA_F7404_US1_GV269_0002.cut.vrt

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