#!/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 1695.0 1068.0 21 56.75 -gcp 9930.0 1032.0 22.5 56.75 -gcp 10008.0 7386.0 22.5 56.11666666666667 -gcp 1656.0 7407.0 21 56.11666666666667 k_001_kta-3-415_32.jpeg k_001_kta-3-415_32.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1651.5 7121.921710161022\n1155.0 7117.519988895059\n1185.0 2739.015824541921\n1687.5 2730.0508051082734\n1695.0 1068.0\n9930.0 1032.0\n10008.0 7386.0\n3243.0 7399.500832870627\n3148.5 7606.500832870627\n3111.0 7696.49333703498\n3087.0 7810.546640755136\n2968.5 8052.16490838423\n2869.5 8109.192393114936\n2743.5 8124.2082176568565\n2136.0 8125.531649083843\n2113.5 8100.017490283176\n1815.75 8098.497084952804\n1823.25 7646.940449750138\n1857.0 7647.067046085509\n1856.25 7405.40255413659\n1656.0 7407.0" | \
  gdaltransform -tps -output_xy k_001_kta-3-415_32.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  k_001_kta-3-415_32.vrt k_001_kta-3-415_32.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_32.cut.vrt k_001_kta-3-415_32.xyz