#!/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.345589018088 56.84531944843 -gcp 15690.0 0.0 25.473901112842 56.890635767704 -gcp 15690.0 18660.0 25.556088481431 56.808824310601 -gcp 0.0 18660.0 25.427776386677 56.763507991327 LVVA_F7404_US1_GV2101_0004.jpg LVVA_F7404_US1_GV2101_0004.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"5085.0 1347.0\n5649.0 1491.0\n6225.0 1605.0\n7386.0 1881.0\n8706.0 2016.0\n9552.0 2190.0\n10467.0 2412.0\n10878.0 2481.0\n11700.0 2739.0\n12358.5 2872.5\n13119.0 2958.0\n13500.0 3066.0\n14220.0 3510.0\n14526.0 4020.0\n14658.0 4380.0\n14916.0 5412.0\n15186.0 6888.0\n15342.0 9144.0\n15054.0 10032.0\n14310.0 11358.0\n13194.0 13518.0\n12840.0 13980.0\n12414.0 14430.0\n11196.0 15648.0\n10518.0 16140.0\n9660.0 16647.0\n8754.0 17085.0\n7668.0 17409.0\n7131.0 17517.0\n5994.0 17643.0\n4990.5 17578.5\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV2101_0004.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV2101_0004.vrt LVVA_F7404_US1_GV2101_0004.cut.vrt

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