#!/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.149276733398 57.607624145573 -gcp 12285.0 0.0 27.162237167358 57.548352693748 -gcp 12285.0 21120.0 26.967830657959 57.534026558765 -gcp 0.0 21120.0 26.954870223999 57.59329801059 LVVA_F7404_US1_GV1797_0002.jpg LVVA_F7404_US1_GV1797_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"871.5 16798.5\n780.0 16488.0\n791.25 16184.25\n807.75 16020.75\n828.75 15993.75\n867.0 15966.0\n1050.0 15258.0\n1030.5 10197.0\n870.0 777.0\n11349.0 612.0\n11440.5 10147.5\n11523.0 17835.0\n11617.5 17806.5\n11676.0 17958.0\n11676.0 18162.0\n11482.5 18303.0\n11427.0 18286.5\n11350.5 18295.5\n11298.0 18306.0\n11248.5 18370.5\n11115.0 18357.0\n11013.0 18388.5\n10936.5 18463.5\n10914.0 18553.5\n10795.5 18550.5\n10591.5 18669.0\n10498.5 18658.5\n10470.0 18756.0\n10372.5 18813.0\n10282.5 18829.5\n10177.5 18775.5\n10047.0 18825.0\n9997.5 18739.5\n9795.0 18736.5\n9718.5 18691.5\n6640.5 18109.5\n5154.0 17820.0\n4146.0 17802.0\n3267.0 18367.5\n3112.5 18019.5\n2971.5 17643.0\n2751.0 17464.5\n2541.0 17359.5\n2278.5 17292.0\n1864.5 17236.5\n1356.0 17196.0\n948.0 17113.5\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1797_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1797_0002.vrt LVVA_F7404_US1_GV1797_0002.cut.vrt

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