#!/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 1392.0 1695.0 27.908883333333332 56.924366666666664 -gcp 21927.0 1686.0 29.11703333333333 56.9402 -gcp 22092.0 16356.0 29.1317 56.465333333333334 -gcp 1488.0 16626.0 27.939783333333335 56.44968333333333 1915_3verst_ru_126k/9-7.jpg 9-7.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1437.75 9162.000832870628\n1392.0 1695.0\n11659.5 1688.9991671293724\n21927.0 1686.0\n22015.5 9022.500832870628\n22092.0 16356.0\n11790.0 16486.49750138812\n1488.0 16626.0" | \
  gdaltransform -tps -output_xy 9-7.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  9-7.vrt 9-7.cut.vrt

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