#!/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 26.893672943115 57.443148124767 -gcp 9621.0 0.0 27.025594711304 57.446057890874 -gcp 9621.0 5847.0 27.027740478516 57.405901165845 -gcp 0.0 5847.0 26.896591186523 57.403080684468 LVVA_F7404_US1_GV437_0002.jpg LVVA_F7404_US1_GV437_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"336.0 348.0\n7542.0 366.0\n7717.5 1281.0\n7704.0 1482.0\n7527.0 1593.0\n7528.5 2052.0\n7428.0 2491.5\n7815.0 3141.0\n7950.0 3561.0\n8180.25 3779.25\n8287.5 3928.5\n8274.0 4047.75\n8244.0 4158.0\n8514.0 4189.5\n8605.5 4240.5\n8643.0 4360.5\n8652.0 4521.0\n8589.0 4617.0\n8421.0 4645.5\n8320.5 4644.0\n8230.5 4591.5\n8119.5 4633.5\n8023.5 4626.0\n7954.5 4659.0\n7716.0 4719.0\n7492.5 4707.0\n7453.5 4794.0\n7263.0 4806.0\n7186.5 4701.0\n6774.0 4632.0\n6613.5 4528.5\n6519.0 4336.5\n6442.5 4171.5\n6304.5 4057.5\n6061.5 3969.0\n5767.5 3775.5\n5653.5 3771.0\n5445.0 3534.0\n4986.0 3339.0\n3885.0 3516.0\n492.0 3573.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV437_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV437_0002.vrt LVVA_F7404_US1_GV437_0002.cut.vrt

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