#!/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.86306344 57.25647594 -gcp 22464.0 0.0 26.87422306 57.34332669 -gcp 22464.0 21897.0 27.03146189 57.33831238 -gcp 0.0 21897.0 27.02262658 57.25108517 GX8464/GX_8464_SK_frame_249.jpg GX_8464_SK_frame_249.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "2925.0 20355.0\n2910.0 20037.0\n2724.0 19704.0\n2451.0 19473.0\n2079.0 19350.0\n792.0 19359.0\n744.0 11856.0\n1005.0 11622.0\n1008.0 10224.0\n729.0 9942.0\n675.0 2172.0\n2628.0 2142.0\n2625.0 384.0\n10335.0 330.0\n10605.0 591.0\n12033.0 579.0\n12330.0 279.0\n17025.0 270.0\n19479.0 231.0\n21873.0 198.0\n21975.0 309.0\n22044.0 9726.0\n21798.0 10026.0\n21804.0 11436.0\n22062.0 11739.0\n22146.0 19284.0\n21132.0 19302.0\n20718.0 19413.0\n20373.0 19692.0\n20184.0 20217.0\n20196.0 21270.0\n12381.0 21360.0\n12201.0 21093.0\n10788.0 21111.0\n10563.0 21387.0\n6999.0 21405.0\n5187.0 21408.0\n2940.0 21417.0" | \
  gdaltransform -tps -output_xy GX_8464_SK_frame_249.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_8464_SK_frame_249.vrt GX_8464_SK_frame_249.cut.vrt

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