#!/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 1374.0 1650.0 30.33265 58.370583333333336 -gcp 22032.0 1533.0 31.583716666666668 58.364133333333335 -gcp 22047.0 16392.0 31.56735 57.8898 -gcp 1398.0 16278.0 30.332516666666667 57.8954 1915_3verst_ru_126k/6-9.jpg 6-9.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1387.5012701100761 8979.008328706273\n1374.0 1650.0\n11706.002540220152 1612.511660188784\n22032.0 1533.0\n22037.99872988992 8973.005830094391\n22047.0 16392.0\n11720.998729889923 16321.492504164353\n1398.0 16278.0" | \
  gdaltransform -tps -output_xy 6-9.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  6-9.vrt 6-9.cut.vrt

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