#!/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 1819.6861736334404 531.1907273736813 47.99419535 72.0003961 -gcp 16443.242604501607 733.0494169905608 59.99761516 72.00055957 -gcp 17773.974115755627 16555.198084397554 59.99803244 68.00048691 -gcp 30.823794212218672 16306.669072737368 47.99760554 68.00032752 -gcp 9128.25 994.6203498056636 53.99734188 72.00048184 -gcp 8891.625 16861.182259855635 53.99780697 68.0004111 vlasenko/north/r-39_40.jpg r-39_40.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1819.6861736334404 531.1907273736813\n16443.242604501607 733.0494169905608\n17773.974115755627 16555.198084397554\n30.823794212218672 16306.669072737368" | \
  gdaltransform -tps -output_xy r-39_40.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  r-39_40.vrt r-39_40.cut.vrt

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