#!/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 2014.9116965836713 1050.417732073374 24.75 57.25 -gcp 10134.44991314418 1065.366870483602 26.25 57.25 -gcp 10180.887376954255 7425.363535297387 26.266666666666666 56.61666666666667 -gcp 1951.8621887666473 7399.903557531962 24.75 56.61666666666667 k_001_kta-3-415_41.jpeg k_001_kta-3-415_41.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "2014.9116965836713 1050.417732073374\n10134.44991314418 1065.366870483602\n10180.887376954255 7425.363535297387\n1951.8621887666473 7399.903557531962" | \
  gdaltransform -tps -output_xy k_001_kta-3-415_41.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  k_001_kta-3-415_41.vrt k_001_kta-3-415_41.cut.vrt

# Generate tiles
gdal2tiles.py -r $GDAL_SAMPLING_TILE --xyz -z $GDAL_TILE_ZOOMS -x --processes=$GDAL_TILE_PROCESSES k_001_kta-3-415_41.cut.vrt k_001_kta-3-415_41.xyz