#!/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 0.0 0.0 27.782776639127142 57.35304647830861 -gcp 14520.0 0.0 28.124871460641415 57.35481107818937 -gcp 14520.0 16686.0 28.126391241602036 57.14289667960417 -gcp 0.0 16686.0 27.786047570777413 57.141875351775695 274-Simani_1922.jpg 274-Simani_1922.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"1854.0 1568.809064408662\n1791.75 8144.996252082176\n1752.1875 14735.252915047196\n13286.25 14808.00187395891\n13296.0 13598.625832870628\n12421.875 13590.381871182675\n12473.625 1666.4975013881178\" | \
  gdaltransform -tps -output_xy 274-Simani_1922.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  274-Simani_1922.vrt 274-Simani_1922.cut.vrt

# Generate tiles
gdal2tiles.py -r $GDAL_SAMPLING_TILE --xyz -z $GDAL_TILE_ZOOMS -x --processes=$GDAL_TILE_PROCESSES 274-Simani_1922.cut.vrt 274-Simani_1922.xyz