#!/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 3840.0 3594.0 18 56 -gcp 30336.0 3726.0 21 56 -gcp 31023.0 35286.0 21 54 -gcp 3201.0 35181.0 18 54 N-34-A_DANCIG_XII.1944.jpg N-34-A_DANCIG_XII.1944.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "17052.0 3810.086714841578\n3831.0 3587.859922178988\n3201.0 35181.1834352418\n17083.5 35391.08671484158\n31026.0 35288.923290717066\n30337.5 3724.395775430795" | \
  gdaltransform -tps -output_xy N-34-A_DANCIG_XII.1944.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  N-34-A_DANCIG_XII.1944.vrt N-34-A_DANCIG_XII.1944.cut.vrt

# Generate tiles
gdal2tiles.py -r $GDAL_SAMPLING_TILE --xyz -z $GDAL_TILE_ZOOMS -x --processes=$GDAL_TILE_PROCESSES N-34-A_DANCIG_XII.1944.cut.vrt N-34-A_DANCIG_XII.1944.xyz