#!/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.230531692505 56.753240495206 -gcp 10551.0 0.0 25.31644821167 56.757334172435 -gcp 10551.0 8697.0 25.325546264648 56.714917134297 -gcp 0.0 8697.0 25.240831375122 56.710865920852 LVVA_F7404_US1_GV662_0002.jpg LVVA_F7404_US1_GV662_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"1098.0 4140.0\n967.5 3910.5\n813.0 3769.5\n615.0 3631.5\n609.0 528.0\n10014.0 642.0\n9996.0 7728.0\n8628.0 7717.5\n8428.5 7845.0\n7810.5 7839.0\n7657.5 7495.5\n7561.5 7458.0\n7513.5 7395.0\n7422.0 7347.0\n7350.0 7216.5\n7032.0 7188.0\n6780.0 7236.0\n2166.0 5586.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV662_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV662_0002.vrt LVVA_F7404_US1_GV662_0002.cut.vrt

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