#!/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 24.37271833 57.07336736 -gcp 10800.0 0.0 24.38767433 56.98837061 -gcp 10800.0 10800.0 24.22811508 56.98057234 -gcp 0.0 10800.0 24.21378136 57.06575023 GX220/GX_220_SK_frame_558.jpg GX_220_SK_frame_558.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "156.0 5760.0\n276.0 5631.0\n291.0 4935.0\n147.0 4758.0\n117.0 879.0\n1044.0 888.0\n1041.0 300.0\n3318.0 321.0\n3330.0 9.0\n4950.0 9.0\n5037.0 156.0\n5787.0 162.0\n5913.0 21.0\n8835.0 24.0\n8835.0 345.0\n10731.0 396.0\n10720.5 9567.0\n10360.5 9555.0\n10126.5 9598.5\n9942.0 9748.5\n9858.0 10050.0\n9849.0 10533.0\n5937.0 10548.0\n5823.0 10416.0\n5101.5 10407.0\n4953.0 10548.0\n1140.0 10543.5\n1152.0 9978.0\n1108.5 9768.0\n1002.0 9633.0\n846.0 9522.0\n630.0 9487.5\n138.0 9480.0" | \
  gdaltransform -tps -output_xy GX_220_SK_frame_558.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_220_SK_frame_558.vrt GX_220_SK_frame_558.cut.vrt

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