#!/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 26.53244745552639 57.352387065286166 -gcp 14340.0 0.0 26.870850249494847 57.3531367087284 -gcp 14340.0 16344.0 26.870409407115137 57.1454543426198 -gcp 0.0 16344.0 26.534413062946303 57.14501670800011 215-Vecgulbene_1924.jpg 215-Vecgulbene_1924.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"7174.875 1519.1229178234314\n1868.4375 1502.6247917823432\n1830.0 8067.749791782342\n1812.0 14634.0\n7126.875 14647.874167129372\n12475.125 14664.562812326483\n12471.0 1543.311563020544\" | \
  gdaltransform -tps -output_xy 215-Vecgulbene_1924.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  215-Vecgulbene_1924.vrt 215-Vecgulbene_1924.cut.vrt

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