#!/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 1078.9504823151126 634.9272626318713 41.99827759 56.00003046 -gcp 14350.865594855306 667.8931149361466 47.99839502 56.00011147 -gcp 14992.778778135047 16434.32731815658 47.99854191 52.0000371 -gcp 392.0594855305467 16420.848972792894 41.99843522 51.99996004 -gcp 3292.5 716.3503609106052 42.99891802 55.97988962 -gcp 5499.75 764.3470294280955 43.99595678 55.96861469 -gcp 7721.25 793.6249305941144 44.9993181 55.96231484 -gcp 9925.5 783.1457523598001 45.99525101 55.96655554 -gcp 12136.875 744.8114242087729 46.99537811 55.97820135 -gcp 14502.75 4612.250347029428 47.99458931 54.96380226 -gcp 14662.5 8554.98431426985 47.99438471 53.95181038 -gcp 14823.9375 12493.728831204884 47.99496855 52.96447559 -gcp 12559.875 16510.223486951694 46.99668685 51.9811373 -gcp 10127.25 16560.97737368129 45.99646216 51.96832726 -gcp 7694.25 16573.607023875626 44.99715434 51.96473054 -gcp 5260.125 16558.606191004998 43.99786258 51.96780467 -gcp 2828.25 16505.72910882843 42.99937082 51.98002044 -gcp 558.75 12475.031302054416 41.99629768 52.96438855 -gcp 732.0 8528.263742365352 41.99686037 53.95248994 -gcp 905.625 4577.6943364797335 41.9975428 54.96530207 -gcp 7703.25 4737.443711826763 44.99377938 54.92758994 -gcp 10023.0 8679.867018323153 45.99439051 53.91812838 -gcp 7693.875 12628.706968350918 44.99466399 52.92872263 -gcp 5375.625 8669.01075791227 43.99515089 53.9188623 -gcp 7696.5 8689.6047334814 44.99334912 53.91465282 vlasenko/north/n-38.jpg n-38.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1078.9504823151126 634.9272626318713\n14350.865594855306 667.8931149361466\n14992.778778135047 16434.32731815658\n392.0594855305467 16420.848972792894" | \
  gdaltransform -tps -output_xy n-38.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  n-38.vrt n-38.cut.vrt

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