#!/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 4346.25 1202.400749583565 35.98053451 84.00038772 -gcp 14185.5 1186.983273181566 59.99295016 84.00151749 -gcp 17437.5 16725.167406996115 59.99575681 80.00069575 -gcp 1083.75 16751.46404775125 35.99414136 80.00032847 -gcp 5973.375 1488.1047334813993 39.99060722 84.00045733 -gcp 7614.75 1661.472272348695 43.99099258 84.00052413 -gcp 9267.375 1717.1960022209882 47.99142183 84.00058797 -gcp 10918.5 1660.7257079400333 51.99189289 84.00064855 -gcp 12558.0 1476.1232648528594 55.99240346 84.00070556 -gcp 14998.875 5070.869100499723 59.99395409 83.00074333 -gcp 15811.875 8956.942046085509 59.99470579 82.00072769 -gcp 16626.0 12840.242157134926 59.99528994 81.00071183 -gcp 14735.25 17190.107856746254 55.99542716 80.00064309 -gcp 12007.875 17483.34661299278 51.9951198 80.00058663 -gcp 9259.125 17581.206829539144 47.99483622 80.00052663 -gcp 6517.125 17493.82599944475 43.9945778 80.0004634 -gcp 3788.25 17216.306635202665 39.9943458 80.00039724 -gcp 1903.5 12873.13784008884 35.99349678 81.0003435 -gcp 2716.875 8983.53768739589 35.99269027 82.00035842 -gcp 3532.6875 5089.736327040533 35.99165245 83.00037321 -gcp 9268.5 5684.433647973348 47.99264244 83.00057291 -gcp 9268.125 9648.881177123822 47.99355721 82.00055766 -gcp 9267.0 13619.48112159911 47.99426808 81.00054223 -gcp 7068.0 9579.813298167684 43.99323479 82.00049408 -gcp 11465.53125 9576.397383398113 51.99391102 82.00061798 vlasenko/north/u-37_38_39_40.jpg u-37_38_39_40.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "4346.25 1202.400749583565\n14185.5 1186.983273181566\n17437.5 16725.167406996115\n1083.75 16751.46404775125" | \
  gdaltransform -tps -output_xy u-37_38_39_40.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  u-37_38_39_40.vrt u-37_38_39_40.cut.vrt

# Generate tiles
gdal2tiles.py -r $GDAL_SAMPLING_TILE --xyz -z $GDAL_TILE_ZOOMS -x --processes=$GDAL_TILE_PROCESSES u-37_38_39_40.cut.vrt u-37_38_39_40.xyz