#!/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 1333.4999999993888 1343.2946081146888 26.684761047363285 56.4042137774986 -gcp 9088.499999999229 1341.17259588512 28.20877075195313 56.362966816267914 -gcp 8253.0 4740.353529738744 28 56 -gcp 3078.0 5028.530294608116 27 56 k_001_kta-3-414_25.jpeg k_001_kta-3-414_25.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1333.5 1343.2946081156197\n9088.5 1341.1725958866036\n9375.0 1790.8198999444137\n9477.0 1958.788215675375\n9459.0 2156.851584213452\n9288.0 2565.023346303502\n9237.0 2746.4824902723735\n9084.0 2889.0225125069483\n9082.5 7906.88938299055\n1317.0 7923.358532518067" | \
  gdaltransform -tps -output_xy k_001_kta-3-414_25.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  k_001_kta-3-414_25.vrt k_001_kta-3-414_25.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_25.cut.vrt k_001_kta-3-414_25.xyz