#!/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 27.271499633789 57.336991615104 -gcp 14484.0 0.0 27.406425476074 57.331849451367 -gcp 14484.0 11457.0 27.401962280273 57.264193174328 -gcp 0.0 11457.0 27.266607284546 57.269159172617 LVVA_F7404_US1_GV1765_0002.jpg LVVA_F7404_US1_GV1765_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"591.0 6000.0\n2730.0 5961.0\n2745.0 5496.0\n3204.0 4716.0\n4089.0 4212.0\n5019.0 3183.0\n4840.5 2773.5\n4711.5 2481.0\n4849.5 2347.5\n5055.0 2289.0\n5135.25 2242.5\n5207.25 2199.75\n5562.75 1871.25\n5570.25 1823.25\n5547.75 1732.5\n5583.75 1617.0\n5677.5 1449.75\n5810.25 1298.25\n5996.25 1165.5\n6175.5 1104.0\n6368.25 1048.5\n6480.75 1003.5\n6742.5 886.5\n7036.5 710.25\n7168.5 687.75\n7323.0 720.0\n7472.25 664.5\n7627.5 616.5\n9009.0 1716.0\n9142.5 1827.0\n9276.0 1938.0\n9111.0 2046.0\n9105.0 3091.5\n13728.0 3192.0\n13695.0 10242.0\n582.0 10179.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1765_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1765_0002.vrt LVVA_F7404_US1_GV1765_0002.cut.vrt

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