#!/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 1926.0 530.3228761799 59.99695408 76 -gcp 13384.5 537.1911438089951 71.99774808 76.00076003 -gcp 14973.888658500688 16302.103838580722 72 72 -gcp 327.0 16283.343697945584 60 72 -gcp 3990.0 16557.13325930039 62.99776179 72.00059517 -gcp 7655.25 16655.390338700723 65.99791457 72.00062847 -gcp 11326.5 16569.31399222654 68.99807306 72.00065937 -gcp 10525.5 750.0995280399778 68.99753891 76.00073105 -gcp 7655.25 822.8095502498611 65.99733648 76.00069954 -gcp 4783.125224305958 748.8858303891348 62.99714136 76.00066558 vlasenko/north/s-41_42.jpg s-41_42.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1926.0 530.3228761799\n13384.5 537.1911438089951\n14973.888658500688 16302.103838580722\n327.0 16283.343697945584" | \
  gdaltransform -tps -output_xy s-41_42.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  s-41_42.vrt s-41_42.cut.vrt

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