#!/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 2478.75 375.1519988895058 47.99258628 76.00046255 -gcp 13951.5 763.8248195446973 59.99695408 76.00062927 -gcp 14974.5 16530.438922820656 59.99761516 72.00055957 -gcp 339.0 16013.386590782899 47.99709773 72.00039614 -gcp 8205.0 848.3049694614102 53.99660505 76.00055 -gcp 7645.5 16623.929067184898 53.99734188 72.00048184 -gcp 5334.0 679.6049416990561 50.99644427 76.00050726 -gcp 11087.25 883.6365907828983 56.99677515 76.00059071 -gcp 11313.0 16671.240283176016 56.99747506 72.00052176 -gcp 3987.75 16405.655330372017 50.99721599 72.00043993 vlasenko/north/s-39_40.jpg s-39_40.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "2478.75 375.1519988895058\n13951.5 763.8248195446973\n14974.5 16530.438922820656\n339.0 16013.386590782899" | \
  gdaltransform -tps -output_xy s-39_40.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  s-39_40.vrt s-39_40.cut.vrt

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