#!/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 26.026954650879 56.804449885889 -gcp 15945.0 0.0 26.179389953613 56.801160291035 -gcp 15945.0 15945.0 26.169605255127 56.717130949426 -gcp 0.0 15945.0 26.022663116455 56.718072968735 LVVA_F7404_US1_GV1114_0002.jpg LVVA_F7404_US1_GV1114_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"779.25 11551.5\n579.0 11409.75\n580.5 10990.5\n774.0 10873.5\n795.0 5646.0\n789.0 417.0\n14754.0 405.0\n14761.5 3036.0\n15103.5 3172.5\n15129.0 4317.0\n14781.0 5128.5\n14844.0 10293.0\n15201.0 10554.0\n15246.0 11082.0\n15249.0 11496.0\n15225.0 12333.0\n15195.0 12612.0\n14853.0 12813.0\n14892.0 14928.0\n10434.0 14844.0\n5508.0 14850.0\n729.0 14910.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1114_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1114_0002.vrt LVVA_F7404_US1_GV1114_0002.cut.vrt

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