#!/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.88844062 56.84665084 -gcp 21435.0 0.0 26.89728118 56.94279034 -gcp 21435.0 21381.0 27.06301354 56.93789643 -gcp 0.0 21381.0 27.0626273 56.84819948 GX825/GX_825_SK_frame_127.jpg GX_825_SK_frame_127.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "408.0 11790.0\n714.0 11574.0\n732.0 10068.0\n414.0 9714.0\n390.0 2040.0\n2259.0 2010.0\n2283.0 426.0\n9894.0 510.0\n10218.0 630.0\n11532.0 648.0\n12024.0 546.0\n21348.0 462.0\n21393.0 9819.0\n21163.5 10110.0\n21138.0 10782.0\n21169.5 11499.0\n21366.0 11790.0\n21312.0 19518.0\n20790.0 19686.0\n20256.0 20082.0\n19950.0 20592.0\n19788.0 21246.0\n11775.0 21189.0\n11598.0 21006.0\n10189.5 21007.5\n9832.5 21171.0\n5799.0 21192.0\n2868.0 21222.0\n2772.0 20520.0\n2448.0 19884.0\n1986.0 19368.0\n1236.0 18984.0\n354.0 18528.0" | \
  gdaltransform -tps -output_xy GX_825_SK_frame_127.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_825_SK_frame_127.vrt GX_825_SK_frame_127.cut.vrt

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