#!/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 1834.8300521134915 1134.2392996108947 27.25 57.75 -gcp 8135.0489287782275 1138.5575319622012 28.383333333333333 57.75 -gcp 8162.094093804284 7485.2367982212345 28.383333333333333 57.11666666666667 -gcp 1776.3317892298785 7476.328515842134 27.25 57.11666666666667 k_001_kta-3-415_43.jpeg k_001_kta-3-415_43.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1834.8300521134915 1134.2392996108947\n2467.511291256514 1136.2379099499722\n2464.504342790967 1985.2529182879377\n3436.5686160972787 2826.0266814897163\n4671.020845396642 2805.0\n8268.013896931094 5416.522512506948\n8260.557324840764 7485.0\n1776.3317892298785 7476.328515842134" | \
  gdaltransform -tps -output_xy k_001_kta-3-415_43.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  k_001_kta-3-415_43.vrt k_001_kta-3-415_43.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_43.cut.vrt k_001_kta-3-415_43.xyz