#!/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.28013142952008 57.449927380640275 -gcp 14340.0 0.0 27.618576098053232 57.44968031363106 -gcp 14340.0 10731.0 27.617439621958948 57.313402117203 -gcp 0.0 10731.0 27.27965977950548 57.314055174402704 250-Lipna_1923.jpg 250-Lipna_1923.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"7285.5 1289.9791782343143\n1978.125 1308.0\n1983.9375 5257.687916435314\n1996.125 9205.688436979455\n7307.25 9184.122917823432\n12609.0 9169.872709605774\n12592.3125 5229.93812465297\n12573.375 1287.7487506940588\" | \
  gdaltransform -tps -output_xy 250-Lipna_1923.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  250-Lipna_1923.vrt 250-Lipna_1923.cut.vrt

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