#!/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 2029.9369774919614 679.8847862298724 11.99636856 71.99977182 -gcp 16660.95434083601 684.3772903942255 23.99646051 71.99999347 -gcp 18208.479656067044 16416.08902727405 23.9970798 67.99993477 -gcp 465.2392282958199 16422.403942254303 11.99700393 67.99971859 -gcp 9348.75 1038.1873958911715 17.99639478 71.99988335 -gcp 9342.375 16848.908661854526 17.99702557 67.99982737 -gcp 4470.75 887.3649361465851 13.99637288 71.99980908 -gcp 6908.25 1002.1807329261522 15.99638163 71.99984627 -gcp 11797.5 999.3323153803442 19.99641233 71.99992027 -gcp 14237.26286587779 880.7825589136327 21.99643425 71.99995699 -gcp 17049.0 4617.420599666852 23.99664033 70.99997881 -gcp 17433.0 8547.224042198777 23.99680182 69.99996413 -gcp 17823.0 12478.69280955025 23.9969476 68.99994946 -gcp 15261.0 16654.978067740143 21.99705813 67.9998992 -gcp 12301.5 16806.446418656302 19.99704005 67.99986338 -gcp 6379.5 16812.37229317046 15.99701472 67.9997912 -gcp 3418.5 16663.818989450305 13.9970075 67.99975493 -gcp 857.625 12484.188506385342 11.9968683 68.99973181 -gcp 1250.25 8540.710785674626 11.99671874 69.99974508 -gcp 1641.0 4604.977096057745 11.99655305 70.99975842 -gcp 9339.75 8947.559550249862 17.99674243 69.99985531 vlasenko/north/r-33_34.jpg r-33_34.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "2029.9369774919614 679.8847862298724\n16660.95434083601 684.3772903942255\n18208.479656067044 16416.08902727405\n465.2392282958199 16422.403942254303" | \
  gdaltransform -tps -output_xy r-33_34.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  r-33_34.vrt r-33_34.cut.vrt

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