#!/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 275.25 316.7456968350916 71.99944818 7.99930292 -gcp 15925.5 288.9219877845641 77.99955509 7.99931017 -gcp 16007.25 10797.953636868408 77.99955832 3.99922814 -gcp 246.375 10822.226193781233 71.99945219 3.99922451 -gcp 8099.25 334.5849528039978 74.99950095 7.99930674 -gcp 8124.0 10822.62722098834 74.99950458 3.99922642 -gcp 13380.75 10807.696557468073 76.99954026 3.99922759 -gcp 13316.25 302.2729039422543 76.9995369 7.99930907 -gcp 2883.0 326.6660188784009 72.99946561 7.99930423 -gcp 2872.5 10822.262423653527 72.99946949 3.99922517 -gcp 261.75 5571.832245974459 71.99945052 5.99926327 -gcp 8110.125 5574.434342032205 74.99950307 5.99926614 -gcp 15966.375 5541.037583287063 77.99955698 5.99926871 vlasenko/north/b-43.jpg b-43.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "275.25 316.7456968350916\n8103.0 331.49728424864213\n15925.5 288.9219877845641\n16007.25 10797.953636868408\n8125.125 10822.132242003621\n246.375 10822.226193781233" | \
  gdaltransform -tps -output_xy b-43.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  b-43.vrt b-43.cut.vrt

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