#!/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 350.25 304.6365907828984 53.99914314 15.99940908 -gcp 15570.75 286.67740144364245 59.99923123 15.9994318 -gcp 15717.75 10801.694475291504 59.99924442 11.99935639 -gcp 247.5 10803.1882287618 53.99915784 11.99933925 -gcp 2885.25 335.67538867295946 54.99915719 15.99941304 -gcp 7957.6875 354.0389367018323 56.99918607 15.99942075 -gcp 13031.4375 316.90935591338143 58.99921594 15.99942819 -gcp 13144.5 10821.24309411438 58.99922939 11.99935366 -gcp 7984.875 10838.676950305386 56.99920003 11.99934805 -gcp 2826.9375 10825.72862298723 54.99917165 11.99934224 -gcp 298.875 5554.939339255969 53.99915107 13.99937381 -gcp 15646.125 5550.027588839534 59.99923835 13.99939375 -gcp 7971.75 5595.843073292615 56.9991936 13.99938405 vlasenko/north/d-40.jpg d-40.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "350.25 304.6365907828984\n7961.25 356.2837312604109\n15570.75 286.67740144364245\n15717.75 10801.694475291504\n7986.75 10838.64498889506\n247.5 10803.1882287618" | \
  gdaltransform -tps -output_xy d-40.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  d-40.vrt d-40.cut.vrt

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