#!/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.284519195557 56.824651065123 -gcp 15924.0 0.0 25.417642593384 56.867976364538 -gcp 15924.0 18567.0 25.494203567505 56.784660969635 -gcp 0.0 18567.0 25.361080169678 56.74133567022 LVVA_F7404_US1_GV2101_0003.jpg LVVA_F7404_US1_GV2101_0003.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"7.5 5326.5\n909.0 5043.0\n2025.0 4377.0\n4455.0 3414.0\n4494.0 3552.0\n4656.0 3438.0\n4698.0 2934.0\n5181.0 2208.0\n5406.0 1830.0\n5562.0 1425.0\n6438.0 312.0\n6894.0 198.0\n7902.0 453.0\n8976.0 660.0\n11334.0 912.0\n12495.0 1035.0\n14394.0 1389.0\n15153.0 1503.0\n15912.0 1539.0\n15921.0 16885.5\n15130.5 17100.0\n14373.0 17203.5\n13950.0 17260.5\n12508.5 17203.5\n12507.0 17616.0\n11313.0 17787.0\n9645.0 17745.0\n8754.0 17652.0\n8016.0 17490.0\n6339.0 16920.0\n5196.0 16566.0\n4212.0 16365.0\n3375.0 16035.0\n3033.0 15708.0\n2658.0 15486.0\n1269.0 15036.0\n3.0 14475.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV2101_0003.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV2101_0003.vrt LVVA_F7404_US1_GV2101_0003.cut.vrt

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