#!/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 1320.0 1187.9291101055805 23.683090209960938 56.737664083138206 -gcp 9118.5 1165.7149321266968 25.208816528320316 56.71731992351671 -gcp 9168.0 7773.592760180995 25.19302368164063 56.02102914268314 -gcp 1389.0 7770.316742081448 23.692016601562504 55.99761259485344 k_001_kta-3-414_18.jpeg k_001_kta-3-414_18.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1317.0 1185.7488687782804\n9118.5 1165.500754147813\n9168.0 7771.4992458521865\n1387.5 7771.502262443439" | \
  gdaltransform -tps -output_xy k_001_kta-3-414_18.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  k_001_kta-3-414_18.vrt k_001_kta-3-414_18.cut.vrt

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