#!/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.860167271481 57.396798628257 -gcp 7104.0 0.0 24.964676813903 57.406226582831 -gcp 7104.0 5853.0 24.976683977936 57.361955590237 -gcp 0.0 5853.0 24.872174435514 57.352527635663 LVVA_F7404_US1_GV257_0002.jpg LVVA_F7404_US1_GV257_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"1834.5 3954.0\n1882.5 3883.5\n1779.0 3774.0\n1443.0 3360.0\n441.0 3351.0\n420.0 384.0\n6756.0 360.0\n6735.0 2068.5\n6018.0 2067.0\n6045.0 2170.5\n6142.5 2248.5\n6162.0 2388.0\n6237.0 2509.5\n6240.0 2653.5\n6301.5 2719.5\n6318.0 2764.5\n6357.0 2791.5\n6396.0 2871.0\n6450.0 2940.0\n6004.5 3603.0\n6412.5 4229.25\n6301.5 4534.5\n6066.75 4663.5\n6012.0 4722.75\n5798.25 4805.25\n5716.5 4731.0\n5671.5 4738.5\n5583.0 4698.0\n2577.75 4496.25\n2439.75 4297.5\n2377.5 4270.5\n2174.25 4348.5\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV257_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV257_0002.vrt LVVA_F7404_US1_GV257_0002.cut.vrt

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