#!/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 288.75 292.7111327040533 156.00094216 11.99931398 -gcp 15786.0 271.67573570238756 162.00100606 11.99929308 -gcp 15901.5 10782.9473903387 162.00099383 7.9992416 -gcp 225.0 10800.9473903387 156.0009307 7.99925559 -gcp 2870.25 314.4111604664075 157.00095355 11.99931059 -gcp 8032.5 325.9323639644642 159.00097545 11.9993037 -gcp 13202.4375 296.1071626873959 161.00099616 11.99929666 -gcp 13281.75 10801.022140477513 161.00098405 7.999244 -gcp 8055.0 10821.593576485286 159.00096359 7.99924871 -gcp 2832.5625 10812.117747084952 157.00094195 7.99925333 -gcp 259.5 5542.9382981676845 156.00093583 9.99928436 -gcp 8041.125 5569.600777345919 159.00096889 9.99927578 -gcp 15842.25 5527.534772348695 162.0009993 9.9992669 vlasenko/north/c-57.jpg c-57.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "288.75 292.7111327040533\n8034.75 324.7737368128818\n15786.0 271.67573570238756\n15901.5 10782.9473903387\n8054.625 10819.51519988895\n225.0 10800.9473903387" | \
  gdaltransform -tps -output_xy c-57.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  c-57.vrt c-57.cut.vrt

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