#!/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 2112.75 450.87159911160467 47.995211 68.00032748 -gcp 19883.25 792.9457245974459 59.9980324 68.00097383 -gcp 21090.0 16537.697806774016 59.99831835 64.00041159 -gcp 296.853027648795 16127.897870318418 47.99795348 64.000257 -gcp 10991.25 1032.8270405330372 53.99780697 68.0004111 -gcp 10678.875 16811.31621321488 53.99812565 64.00033807 vlasenko/north/q-39_40.jpg q-39_40.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "2112.75 450.87159911160467\n19883.25 792.9457245974459\n21090.0 16537.697806774016\n296.853027648795 16127.897870318418" | \
  gdaltransform -tps -output_xy q-39_40.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  q-39_40.vrt q-39_40.cut.vrt

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