#!/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.056892269556 57.419792991232 -gcp 10632.0 0.0 25.210198355432 57.43107188738 -gcp 10632.0 12993.0 25.236673353842 57.334025962258 -gcp 0.0 12993.0 25.083029012157976 57.32270025494399 LVVA_F7404_US1_GV223_0003.jpg LVVA_F7404_US1_GV223_0003.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "3077.25 7991.25\n2889.0 7567.5\n864.0 7614.75\n837.0 2897.25\n1941.0 2875.5\n2007.75 2829.0\n3036.0 2531.25\n3039.0 0.0\n9327.0 0.0\n9405.0 12321.0\n8148.0 12321.0\n8049.0 12498.0\n7240.5 12498.0\n7147.5 12307.5\n3144.0 12291.0" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV223_0003.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  LVVA_F7404_US1_GV223_0003.vrt LVVA_F7404_US1_GV223_0003.cut.vrt

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