#!/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.761810302734 57.868405666727 -gcp 15945.0 0.0 24.908580780029 57.87999899848 -gcp 15945.0 21597.0 24.952783584595 57.765316707384 -gcp 0.0 21597.0 24.797172546387 57.751442372569 LVVA_F7404_US1_GV1939_0002.jpg LVVA_F7404_US1_GV1939_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"5496.0 10152.0\n5476.5 10086.0\n5439.0 10023.0\n5335.5 9904.5\n5073.0 9756.0\n1218.0 9780.0\n1233.0 6613.5\n1366.5 6514.5\n1438.5 6447.0\n1482.0 6404.25\n1515.0 6357.75\n1569.75 6210.0\n1621.5 5988.0\n1647.0 5844.0\n1639.5 5722.5\n1612.5 5671.5\n1579.5 5625.0\n1555.5 5608.5\n1515.0 5599.5\n1473.75 5609.25\n1440.75 5624.25\n1395.0 5675.25\n1356.0 5742.0\n1302.0 5768.25\n1248.0 5761.5\n1205.25 5734.5\n1161.0 5701.5\n1138.5 5648.25\n1119.75 5456.25\n1108.5 5355.0\n1083.0 5274.75\n1052.25 5168.25\n1030.5 5045.25\n1018.5 4937.25\n1099.5 4793.25\n1200.0 4722.75\n1268.25 4715.25\n1423.5 4752.0\n1509.75 4736.25\n1599.0 4734.75\n1680.75 4758.75\n1771.5 4818.0\n1835.25 4819.5\n1915.5 4762.5\n1980.0 4665.75\n1995.75 4551.75\n1986.75 4344.0\n2019.75 4151.25\n2070.0 4051.5\n2136.0 3961.5\n2235.0 3921.0\n2304.75 3879.0\n2466.75 3752.25\n2550.75 3638.25\n2649.75 3454.5\n2803.5 3543.0\n2799.0 2610.0\n3636.0 2604.0\n5469.0 2571.0\n5476.5 1545.0\n14970.0 1500.0\n14970.0 21396.0\n5694.0 21366.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1939_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1939_0002.vrt LVVA_F7404_US1_GV1939_0002.cut.vrt

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