#!/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 1885.5 546.83745141588 71.99774808 76 -gcp 13359.0 504.9140755136035 83.99864051 76.00084917 -gcp 15016.5 16239.936979455857 84 72 -gcp 352.5 16294.65574680733 72 72 -gcp 4016.25 16551.12576346474 74.99840544 72.00071364 -gcp 7684.5 16633.64908384231 77.99857841 72.00073687 -gcp 11356.875 16531.95787062743 80.99875529 72.0007574 -gcp 10505.25 724.6253470294281 80.99841025 76.00083104 -gcp 7629.75 806.344947251527 77.99818435 76.00081009 -gcp 4754.625 748.2016240977234 74.99796342 76.00078641 vlasenko/north/s-43_44.jpg s-43_44.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1885.5 546.83745141588\n13359.0 504.9140755136035\n15016.5 16239.936979455857\n352.5 16294.65574680733" | \
  gdaltransform -tps -output_xy s-43_44.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  s-43_44.vrt s-43_44.cut.vrt

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