#!/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.678812026978 57.535086240698 -gcp 8823.0 0.0 24.765844345093 57.539877461406 -gcp 8823.0 11886.0 24.775285720825 57.476250363545 -gcp 0.0 11886.0 24.687824249268 57.471958459334 LVVA_F7404_US1_GV1888_0002.jpg LVVA_F7404_US1_GV1888_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"433.5 399.0\n7686.0 555.0\n7666.5 4887.0\n8130.0 6921.0\n8526.0 10350.0\n8304.0 10824.0\n7579.5 10917.0\n7047.0 11025.0\n6094.5 11058.0\n4726.5 10467.0\n3634.5 9847.5\n3504.0 9681.0\n3153.0 9459.0\n2670.0 9306.0\n288.0 9666.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1888_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1888_0002.vrt LVVA_F7404_US1_GV1888_0002.cut.vrt

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