#!/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 1866.0 822.1528317601333 23.9970798 67.99993477 -gcp 19716.0 856.7065519156025 35.9972833 68.00014075 -gcp 21219.75 16595.405330372017 35.99767807 64.00007586 -gcp 339.0 16541.444058856192 23.99750414 63.99987608 -gcp 10795.5 1261.5991116046641 29.99716602 68.00003961 -gcp 10772.25 17063.4594669628 29.99757783 63.99997777 vlasenko/north/q-35_36.jpg q-35_36.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1866.0 822.1528317601333\n19716.0 856.7065519156025\n21219.75 16595.405330372017\n339.0 16541.444058856192" | \
  gdaltransform -tps -output_xy q-35_36.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  q-35_36.vrt q-35_36.cut.vrt

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