#!/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 1774.8646302250804 505.74902831760136 23.99646051 71.99999347 -gcp 16426.749839228298 390.3106607440311 35.99670716 72.00020465 -gcp 18082.96495176849 16110.37312604109 35.9972833 68.00014075 -gcp 323.45160771704184 16243.730705163798 23.9970798 67.99993477 -gcp 9104.25 819.9286507495837 29.99656501 72.00010096 -gcp 9199.5 16609.997431982232 29.99716602 68.00003961 -gcp 6235.5 16582.7094669628 27.99713378 68.00000502 -gcp 3276.0 16452.14450305386 25.99710503 67.99997006 -gcp 12173.25 16536.119933370348 31.99720172 68.00007379 -gcp 15125.25 16369.954747362575 33.99724083 68.00010752 -gcp 13987.5 614.3749305941144 33.99665568 72.00017058 -gcp 11554.5 758.4065796779566 31.99660828 72.00013601 -gcp 6651.75 799.6103553581344 27.99652593 72.0000655 -gcp 4213.5 695.6724736257635 25.99649108 72.00002965 -gcp 1046.625 8373.450999444753 23.99680182 69.99996413 -gcp 9153.0 8713.577665186009 29.99689625 70.00007037 -gcp 17255.625 8254.175805108272 35.99702469 70.00017284 -gcp 1410.0 4446.938228761799 23.99664033 70.99997881 -gcp 16841.25 4332.52137701277 35.99687445 71.00018878 -gcp 17668.5 12179.724319822322 35.99716031 69.00015683 -gcp 685.5 12302.368684064408 23.9969476 68.99994946 vlasenko/north/r-35_36.jpg r-35_36.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1774.8646302250804 505.74902831760136\n16426.749839228298 390.3106607440311\n18082.96495176849 16110.37312604109\n323.45160771704184 16243.730705163798" | \
  gdaltransform -tps -output_xy r-35_36.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  r-35_36.vrt r-35_36.cut.vrt

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