#!/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.688940048218 57.697221926924 -gcp 9864.0 0.0 24.83304977417 57.70598153103 -gcp 9864.0 10638.0 24.84712600708 57.631756019121 -gcp 0.0 10638.0 24.700613021851 57.622335004396 LVVA_F7404_US1_GV323_0002.jpg LVVA_F7404_US1_GV323_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"477.0 6951.0\n432.0 405.0\n6783.0 456.0\n6742.5 5521.5\n7059.0 5514.0\n7089.0 9078.0\n6714.0 9082.5\n6453.0 8964.0\n4947.0 9250.5\n3813.0 9283.5\n1699.125 8798.625\n964.5 8451.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV323_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV323_0002.vrt LVVA_F7404_US1_GV323_0002.cut.vrt

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