#!/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 25.981645388306 57.060179822088 -gcp 15069.0 0.0 26.126346588135 57.076574722762 -gcp 15069.0 15936.0 26.154155731201 56.994191940781 -gcp 0.0 15936.0 26.012620925903 56.978385719862 LVVA_F7404_US1_GV945_0002.jpg LVVA_F7404_US1_GV945_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"2412.0 12874.5\n2352.0 12807.0\n2292.0 12687.0\n2269.5 12565.5\n2236.5 12466.5\n2100.0 12076.5\n1983.0 11880.0\n1942.5 11725.5\n1875.0 11568.0\n1801.5 11467.5\n1746.0 11343.0\n1516.5 10885.5\n1477.5 10854.0\n1452.0 10813.5\n1452.0 10747.5\n1231.5 10188.0\n1201.5 10024.5\n1222.5 9891.0\n1218.0 9796.5\n1242.0 9715.5\n1240.5 9585.0\n1141.5 9363.0\n1042.5 9141.0\n918.0 9067.5\n741.0 8986.5\n701.25 8931.0\n696.75 8856.0\n681.0 8840.25\n690.75 8770.5\n741.0 8610.0\n708.0 681.0\n14185.5 726.0\n14306.25 15023.25\n3636.0 14979.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV945_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV945_0002.vrt LVVA_F7404_US1_GV945_0002.cut.vrt

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