#!/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 0.0 0.0 27.03217149 57.43617298 -gcp 21315.0 0.0 27.05670297 57.4281944 -gcp 21315.0 21306.0 27.04254091 57.41518468 -gcp 0.0 21306.0 27.01759636 57.42277328 GX3778/GX_3778_SD_frame_183.jpg GX_3778_SD_frame_183.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "39.0 11658.0\n201.0 11463.0\n303.0 10719.0\n249.0 10062.0\n39.0 9816.0\n49.5 6049.5\n102.0 2289.0\n1272.0 2298.0\n1854.0 2022.0\n2142.0 1566.0\n2220.0 942.0\n2220.0 150.0\n9567.0 195.0\n9846.0 528.0\n11280.0 525.0\n11538.0 210.0\n16473.0 234.0\n19350.0 243.0\n19357.5 2254.5\n20025.0 2334.0\n21036.0 2304.0\n20946.0 9597.0\n20703.0 10170.0\n20694.0 11616.0\n20931.0 12012.0\n20886.0 17202.0\n20880.0 21276.0\n11332.5 21298.5\n11152.5 21079.5\n9768.0 21061.5\n9531.0 21292.5\n69.0 21288.0" | \
  gdaltransform -tps -output_xy GX_3778_SD_frame_183.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_3778_SD_frame_183.vrt GX_3778_SD_frame_183.cut.vrt

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