#!/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.72516584 56.9800222 -gcp 21123.0 0.0 26.89006235 56.99281202 -gcp 21123.0 21498.0 26.9091198 56.89813167 -gcp 0.0 21498.0 26.74306397 56.89050611 GX825/GX_825_SK_frame_105.jpg GX_825_SK_frame_105.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "0.0 2463.0\n939.0 2463.0\n1323.0 2274.0\n1674.0 2025.0\n1866.0 1650.0\n1938.0 1059.0\n1896.0 288.0\n9399.0 366.0\n9651.0 699.0\n11058.0 693.0\n11263.5 402.0\n15180.0 393.0\n19173.0 426.0\n19194.0 2316.0\n20934.0 2364.0\n20904.0 6547.5\n20865.0 10056.0\n20740.5 10201.5\n20574.0 10314.0\n20572.5 11737.5\n20727.0 11833.5\n20847.0 11976.0\n20856.0 13578.0\n20833.5 17559.0\n20797.5 19569.0\n20781.0 21495.0\n11050.5 21495.0\n10911.0 21235.5\n9519.0 21238.5\n9372.0 21495.0\n1734.0 21495.0\n762.0 20856.0\n0.0 19800.0" | \
  gdaltransform -tps -output_xy GX_825_SK_frame_105.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_825_SK_frame_105.vrt GX_825_SK_frame_105.cut.vrt

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