#!/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.78994561 56.98109822 -gcp 21252.0 0.0 26.95491062 56.99364808 -gcp 21252.0 21498.0 26.97496528 56.89804663 -gcp 0.0 21498.0 26.80566582 56.89201604 GX825/GX_825_SK_frame_106.jpg GX_825_SK_frame_106.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1332.0 20898.0\n672.0 20226.0\n0.0 19386.0\n0.0 11691.0\n156.0 11559.0\n150.0 10146.0\n0.0 10035.0\n0.0 2394.0\n993.0 2445.0\n1473.0 2250.0\n1902.0 1821.0\n2043.0 1210.5\n2049.0 267.0\n6372.0 261.0\n9516.0 288.0\n9654.0 381.0\n9783.0 612.0\n11169.0 618.0\n11337.0 390.0\n11523.0 312.0\n19329.0 342.0\n19353.0 2247.0\n21102.0 2274.0\n21045.0 9999.0\n20718.0 10230.0\n20724.0 11667.0\n20907.0 11766.0\n21030.0 11934.0\n21015.0 16761.0\n20970.0 19149.0\n20955.0 21495.0\n11307.0 21478.5\n11199.0 21339.0\n11127.0 21180.0\n10395.0 21153.0\n9684.0 21171.0\n9579.0 21345.0\n9459.0 21495.0\n1659.0 21453.0" | \
  gdaltransform -tps -output_xy GX_825_SK_frame_106.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_825_SK_frame_106.vrt GX_825_SK_frame_106.cut.vrt

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