#!/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 1816.581185087351 5512.014467428998 25.20179987 57.17847074 -gcp 19702.32152821848 6219.153766969604 25.20316243 57.10749705 -gcp 19280.326064685793 20120.96784650944 25.10309458 57.10532951 -gcp 3312.4185963528084 20243.001765196022 25.09346008 57.16942115 GX10514/GX_10514_SK_frame_436.jpg GX_10514_SK_frame_436.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1248.0 3888.0493573264794\n20688.0 3647.7778920308483\n20736.0 20688.172750642672\n1008.0 20929.33264781491" | \
  gdaltransform -tps -output_xy GX_10514_SK_frame_436.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_10514_SK_frame_436.vrt GX_10514_SK_frame_436.cut.vrt

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