#!/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 24.152498245239 56.88265709052 -gcp 10596.0 0.0 24.255065917969 56.88706487309 -gcp 10596.0 8856.0 24.266910552979 56.841672669977 -gcp 0.0 8856.0 24.164342880249 56.837264887407 LVVA_F7404_US1_GV82_0002.jpg LVVA_F7404_US1_GV82_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"4803.0 6066.0\n4500.0 5556.0\n4073.25 4908.75\n3965.25 4796.25\n3867.0 4666.5\n3622.5 4281.0\n3435.0 4002.0\n3337.5 3721.5\n3216.0 3117.0\n3219.0 2913.0\n3336.0 2406.0\n3417.0 2169.0\n3837.0 1602.0\n4296.0 1218.0\n5286.0 639.0\n5475.0 516.75\n5501.25 394.5\n5580.0 369.75\n6085.5 328.5\n6276.0 342.0\n6394.5 378.75\n6417.0 598.5\n6444.75 750.0\n6516.0 948.0\n6664.5 1185.0\n6763.5 1332.0\n6870.0 1453.5\n7281.0 1728.0\n7686.0 2082.0\n7848.0 2190.0\n8074.5 2173.5\n8308.5 2178.0\n8424.0 2220.0\n8551.5 2320.5\n8653.5 2434.5\n8854.5 2709.0\n9057.0 2835.0\n9240.0 2988.0\n9397.5 3225.0\n9472.5 3411.0\n9547.5 3492.0\n9615.0 3679.5\n9775.5 3973.5\n9772.5 4122.0\n9832.5 4215.0\n9904.5 4357.5\n10015.5 4470.0\n10011.0 7848.0\n7008.0 7804.5\n6940.5 8049.0\n5946.0 8059.5\n5837.625 7801.5\n5656.5 7519.5\n5460.0 7060.5\n5335.5 6787.5\n5239.5 6633.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV82_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV82_0002.vrt LVVA_F7404_US1_GV82_0002.cut.vrt

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