#!/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 284.25 321.8574403109384 138.00069188 11.99936695 -gcp 15760.5 272.40991116046644 144.00078421 11.99935111 -gcp 15896.25 10779.93614658523 144.00077467 7.99928046 -gcp 233.25 10818.959050527485 138.00068347 7.99929106 -gcp 8021.25 338.3232926152138 141.00073906 11.99935928 -gcp 8058.375 10820.318087173791 141.00073007 7.99928592 -gcp 13287.375 10792.211410327596 143.00076003 7.99928232 -gcp 13185.5625 300.02748473070517 143.00076939 11.99935389 -gcp 2864.8125 337.9097723486952 139.00070782 11.99936445 -gcp 2841.1875 10829.130899500276 139.00069922 7.99928939 -gcp 259.6875 5567.679657134925 138.00068723 9.99932861 -gcp 8041.5 5584.520093003886 141.00073409 9.9993222 -gcp 15827.625 5526.993302332037 144.00077894 9.99931538 vlasenko/north/c-54.jpg c-54.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "284.25 321.8574403109384\n8022.001911639762 339.0\n15760.5 272.40991116046644\n15896.25 10779.93614658523\n8059.5 10820.261451971128\n233.25 10818.959050527485" | \
  gdaltransform -tps -output_xy c-54.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  c-54.vrt c-54.cut.vrt

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