#!/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 27.017612457275 57.569094874486 -gcp 15060.0 0.0 27.15507263666 57.581129897304 -gcp 15060.0 12402.0 27.171678543091 57.519648699513 -gcp 0.0 12402.0 27.03364678086195 57.507358515476184 LVVA_F7404_US1_GV1804_0002.jpg LVVA_F7404_US1_GV1804_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"774.0 411.0\n2157.0 417.0\n2949.0 417.0\n3000.0 255.0\n3283.5 244.5\n3381.0 309.0\n3558.0 328.5\n3630.0 391.5\n3693.0 459.0\n5278.5 2116.5\n5749.5 2647.5\n7749.0 3441.0\n7849.5 3495.0\n7950.0 3522.0\n8262.0 3534.0\n8554.5 3645.0\n8644.5 3531.0\n8718.0 3448.5\n8742.0 3355.5\n8799.0 3211.5\n8872.5 3159.0\n9049.5 3195.0\n9186.0 3267.0\n9252.0 3415.5\n9321.0 3484.5\n9387.0 3550.5\n9450.0 3633.0\n9498.0 3717.0\n9522.0 3822.0\n9520.5 3978.0\n9613.5 4074.0\n9826.5 4273.5\n10053.0 4515.0\n10131.0 4561.5\n10179.0 4644.0\n10260.0 4680.0\n10294.5 4773.0\n10441.5 4842.0\n10464.0 4930.5\n10774.5 4932.0\n10275.0 5358.0\n10089.0 5466.0\n10110.0 11232.0\n729.0 11259.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1804_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1804_0002.vrt LVVA_F7404_US1_GV1804_0002.cut.vrt

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