#!/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 756.75 221.41324264297612 41.99807408 60.00010031 -gcp 8660.25 248.45072182121044 47.99820539 60.00018488 -gcp 9084.0 10737.901582454191 47.99839502 56.00011147 -gcp 249.0 10714.67282065519 41.99827759 56.00003046 -gcp 4707.1875 322.04166088284285 44.99813718 60.00014334 -gcp 4663.125 10811.482266796225 44.99833402 56.00007168 -gcp 2072.25 273.04298642533934 42.99747048 59.98114653 -gcp 3390.0 306.42364253393663 43.99745476 59.96954194 -gcp 6025.5 316.5797511312217 45.99654129 59.96910526 -gcp 7345.5 290.10501508295624 46.99868525 59.98158884 -gcp 8764.5 2869.2209653092004 47.99705892 58.9589686 -gcp 8871.0 5493.096719457014 47.99750063 57.94470954 -gcp 8978.25 8118.128205128205 47.99840937 56.95780896 -gcp 7617.75 10772.54600301659 47.00139295 55.98622408 -gcp 6141.23307709031 10801.882164731807 45.99845333 55.97424731 -gcp 3191.1710341104335 10793.02824242853 43.99640547 55.97462126 -gcp 1721.686279137705 10755.892072351959 42.99871714 55.98659115 -gcp 371.4681395688525 8093.240907970851 41.99517945 56.9579331 -gcp 497.25 5472.780542986426 41.99431646 57.94271569 -gcp 627.0 2847.9445701357467 41.99617169 58.95694512 -gcp 4685.625 5570.407616892911 44.99620361 57.91039139 -gcp 4697.625 2943.4204374057317 44.99724206 58.92473243 -gcp 4674.375 8190.1691176470595 44.99572372 56.9265709 -gcp 6081.09375 5559.5120192307695 45.99591974 57.91615905 -gcp 3290.625 5555.4407993966815 43.99700625 57.91447327 vlasenko/north/o-38.jpg o-38.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "756.75 221.41324264297612\n8660.25 248.45072182121044\n9084.0 10737.901582454191\n249.0 10714.67282065519" | \
  gdaltransform -tps -output_xy o-38.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  o-38.vrt o-38.cut.vrt

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