#!/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 2426.25 344.33620210993894 132.00664923 80.0008099 -gcp 13336.5 303.16449194891726 156.00529092 80.0004967 -gcp 15495.75 10688.702873403665 156.00379804 76.0004331 -gcp 302.625 10751.340574680733 132.00238641 76.0007418 -gcp 7903.5 11493.177401443641 144.00316132 76.00060487 -gcp 4075.125 11308.945377568018 138.00278915 76.0006781 -gcp 5142.375 765.6194475291504 138.00388549 80.0007453 -gcp 7876.5 892.9753609106053 144.00440395 80.00067099 -gcp 10615.875 747.4838978345363 150.00487414 80.00058782 -gcp 11719.125 11285.388881177125 150.00349885 76.0005229 -gcp 1894.5 2948.6872570794003 132.0030255 79.00079326 -gcp 13877.25 2899.240352581899 156.00481514 79.00048103 -gcp 14418.0 5494.952456968351 156.00441913 78.0004652 -gcp 14958.5625 8088.970086063298 156.00408448 77.00044922 -gcp 833.0625 8147.285605219323 132.00256639 77.0007592 -gcp 1363.1249785310338 5548.9153259655795 132.00277667 78.00077635 -gcp 7885.5 3542.3570238756247 144.00400792 79.00065477 -gcp 7894.875 6193.1820516379785 144.0036783 78.00063834 -gcp 7896.75 8845.219530816214 144.00339974 77.0006217 -gcp 11172.1875 6016.520301221543 150.00407102 78.00055572 -gcp 4607.4375 6037.347584675181 138.00324526 78.00071215 vlasenko/north/t-53_54_55_56.jpg t-53_54_55_56.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "2426.25 344.33620210993894\n13336.5 303.16449194891726\n15495.75 10688.702873403665\n302.625 10751.340574680733" | \
  gdaltransform -tps -output_xy t-53_54_55_56.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  t-53_54_55_56.vrt t-53_54_55_56.cut.vrt

# Generate tiles
gdal2tiles.py -r $GDAL_SAMPLING_TILE --xyz -z $GDAL_TILE_ZOOMS -x --processes=$GDAL_TILE_PROCESSES t-53_54_55_56.cut.vrt t-53_54_55_56.xyz