#!/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 939.9679732313575 889.6840644086619 21 58 -gcp 11324.574569789675 975.2098833981122 22.5 58 -gcp 11379.970841300192 14085.13450860633 22.5 57 -gcp 713.491395793499 14001.163659078291 21 57 31.Ventspils.1940.jpg 31.Ventspils.1940.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "594.7485659655831 885.7504164353138\n11324.574569789675 975.2098833981122\n11379.970841300192 14085.13450860633\n356.0224665391969 14001.750416435316" | \
  gdaltransform -tps -output_xy 31.Ventspils.1940.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  31.Ventspils.1940.vrt 31.Ventspils.1940.cut.vrt

# Generate tiles
gdal2tiles.py -r $GDAL_SAMPLING_TILE --xyz -z $GDAL_TILE_ZOOMS -x --processes=$GDAL_TILE_PROCESSES 31.Ventspils.1940.cut.vrt 31.Ventspils.1940.xyz