#!/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.05353546 57.07274918 -gcp 21252.0 0.0 24.2091465 57.07398553 -gcp 21252.0 22335.0 24.21313763 56.98627213 -gcp 0.0 22335.0 24.0561533 56.9839514 GX220/GX_220_SK_frame_555.jpg GX_220_SK_frame_555.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "93.0 12738.0\n294.0 12573.0\n279.0 11115.0\n72.0 10914.0\n174.0 3294.0\n1044.0 3297.0\n1314.0 3267.0\n1569.0 3186.0\n1809.0 3009.0\n1980.0 2784.0\n2118.0 2511.0\n2151.0 2241.0\n2085.0 1182.0\n9642.0 1230.0\n9840.0 1482.0\n11268.0 1488.0\n11586.0 1200.0\n19224.0 1224.0\n19224.0 3156.0\n20292.0 3192.0\n20184.0 8304.0\n20952.0 8556.0\n21036.0 22332.0\n11454.0 22332.0\n11274.0 22233.0\n9885.0 22230.0\n9717.0 22332.0\n1974.0 22332.0\n1992.0 21450.0\n1878.0 21138.0\n1650.0 20886.0\n1320.0 20622.0\n900.0 20601.0\n78.0 20604.0" | \
  gdaltransform -tps -output_xy GX_220_SK_frame_555.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_220_SK_frame_555.vrt GX_220_SK_frame_555.cut.vrt

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