#!/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 4517.188085106196 3217.700944969549 22.59749139 57.75510364 -gcp 15208.440640423109 4580.178595241378 27.04767657 57.4231395 -gcp 14140.687659574374 11508.271817676501 26.5244375 55.88571041 -gcp 707.0290069752757 12154.03420241869 21.12808228 55.71434668 -gcp 7927.5 6324.1794342958165 24.02812958 57.06016507 -gcp 10068.0 9039.123747790218 24.89433289 56.44807163 LATVIJAS_KARTE_1921.jpg LATVIJAS_KARTE_1921.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "0.0 0.0\n18408.0 0.0\n18408.0 12810.0\n0.0 12810.0" | \
  gdaltransform -tps -output_xy LATVIJAS_KARTE_1921.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  LATVIJAS_KARTE_1921.vrt LATVIJAS_KARTE_1921.cut.vrt

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