#!/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 25.53371617 56.21629257 -gcp 21435.0 0.0 25.57070171 56.23733949 -gcp 21435.0 21186.0 25.60789533 56.21848289 -gcp 0.0 21186.0 25.57289462 56.1949598 GX4599/GX_4599_SD_frame_21.jpg GX_4599_SD_frame_21.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "2730.0 20412.0\n2668.5 19975.5\n2434.5 19653.0\n2149.5 19419.0\n1818.0 19311.0\n579.0 19296.0\n534.0 15654.0\n508.5 11955.0\n615.0 11709.0\n846.0 11529.0\n861.0 11109.0\n859.5 10557.0\n837.0 10140.0\n604.5 9979.5\n487.5 9781.5\n468.0 8575.5\n411.0 1983.0\n2334.0 1944.0\n2334.0 279.0\n10077.0 165.0\n10368.0 450.0\n11769.0 432.0\n11961.0 165.0\n16587.0 141.0\n21420.0 114.0\n21372.0 10605.0\n21420.0 21147.0\n11973.0 21165.0\n11817.0 21057.0\n11497.5 21010.5\n10861.5 21015.0\n10474.5 21060.0\n10344.0 21174.0\n8847.0 21175.5\n7704.0 21144.0\n6318.0 21162.0\n4896.0 21124.5\n3796.5 21142.5\n2775.0 21142.5" | \
  gdaltransform -tps -output_xy GX_4599_SD_frame_21.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_4599_SD_frame_21.vrt GX_4599_SD_frame_21.cut.vrt

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