#!/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.80231324 57.25824479 -gcp 22077.0 0.0 26.81128254 57.3449725 -gcp 22077.0 21897.0 26.96973333 57.33899624 -gcp 0.0 21897.0 26.96093568 57.25455497 GX8464/GX_8464_SK_frame_250.jpg GX_8464_SK_frame_250.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "537.0 12039.0\n789.0 11817.0\n783.0 10386.0\n507.0 10065.0\n450.0 2331.0\n1476.0 2368.5\n2379.0 2313.0\n2436.0 561.0\n10164.0 501.0\n10377.0 759.0\n11772.0 741.0\n12009.0 474.0\n21549.0 375.0\n21738.0 498.0\n21822.0 9990.0\n21555.0 10209.0\n21567.0 11586.0\n21852.0 11847.0\n21930.0 19528.5\n20844.0 19512.0\n20571.0 19609.5\n20325.0 19756.5\n20122.5 19992.0\n20019.0 20221.5\n20019.0 21474.0\n16326.0 21513.0\n12921.0 21552.0\n12165.0 21579.0\n11973.0 21300.0\n10572.0 21297.0\n10321.5 21588.75\n5556.0 21612.0\n2730.0 21630.0\n2721.0 20508.0\n2703.0 20292.0\n2652.0 20112.0\n2427.0 19815.0\n2121.0 19614.0\n1938.0 19563.0\n1737.0 19536.0\n559.5 19594.5" | \
  gdaltransform -tps -output_xy GX_8464_SK_frame_250.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_8464_SK_frame_250.vrt GX_8464_SK_frame_250.cut.vrt

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