#!/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 26.340923309326175 57.52805951372316 -gcp 10053.0 0.0 26.44619464874268 57.54084484300906 -gcp 10053.0 8463.0 26.466407775878906 57.49830187451385 -gcp 0.0 8463.0 26.35950565338135 57.483448554692615 LVVA_F7404_US1_GV2074_0002.jpg LVVA_F7404_US1_GV2074_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"348.0 219.0\n9711.0 270.0\n9722.25 7370.25\n3379.5 7366.5\n3369.0 7594.5\n3315.0 7681.5\n3231.0 7728.0\n3052.5 7702.5\n2860.5 7659.0\n2680.5 7446.0\n2644.5 7462.5\n2619.0 7485.0\n2520.0 7422.0\n2488.5 7356.0\n312.0 7362.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV2074_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV2074_0002.vrt LVVA_F7404_US1_GV2074_0002.cut.vrt

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