#!/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.779346466064453 57.39905765499924 -gcp 8754.0 0.0 26.790590286254883 57.353849156116205 -gcp 8754.0 18831.0 26.59850120544434 57.34435616327117 -gcp 0.0 18831.0 26.59047603607178 57.39135717729088 LVVA_F7404_US1_GV1807_0002.jpg LVVA_F7404_US1_GV1807_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"1058.25 16061.25\n855.75 15789.0\n815.25 15696.0\n801.0 15610.5\n793.5 15512.25\n834.0 15383.25\n985.5 15012.0\n1062.75 14760.75\n1054.5 13803.75\n1044.0 13764.75\n1033.5 13303.5\n1581.0 12631.5\n1848.0 12417.0\n1878.0 12064.5\n3421.5 10330.5\n3459.0 9276.0\n3477.0 8670.0\n3483.0 8448.0\n3426.0 8206.5\n3427.5 7897.5\n3273.0 7911.0\n3061.5 7864.5\n2928.0 7936.5\n2776.5 7947.0\n2671.5 8007.0\n2416.5 7744.5\n2502.0 7650.0\n2431.5 7350.0\n1225.5 6399.0\n1231.5 6181.5\n1125.0 5983.5\n1143.0 669.0\n7929.0 720.0\n7896.0 4767.0\n7851.0 8805.0\n8322.0 9546.0\n8331.0 9927.0\n7794.0 11349.0\n7784.25 13760.25\n7797.0 13803.75\n7836.0 17841.0\n1074.0 17836.5\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1807_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1807_0002.vrt LVVA_F7404_US1_GV1807_0002.cut.vrt

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