#!/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 21.279691457748 56.644282713813 -gcp 12192.0 0.0 21.28902554512 56.613417843157 -gcp 12192.0 13572.0 21.226626634598 56.607348121806 -gcp 0.0 13572.0 21.217485666275 56.638442184245 Aizpute-1988-991/Aizpute-1988-991-95.jpg Aizpute-1988-991-95.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "2252.25 921.75\n823.5 10731.0\n10585.5 12169.5\n12032.25 2346.75" | \
  gdaltransform -tps -output_xy Aizpute-1988-991-95.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  Aizpute-1988-991-95.vrt Aizpute-1988-991-95.cut.vrt

# Generate tiles
gdal2tiles.py -r $GDAL_SAMPLING_TILE --xyz -z $GDAL_TILE_ZOOMS -x --processes=$GDAL_TILE_PROCESSES Aizpute-1988-991-95.cut.vrt Aizpute-1988-991-95.xyz