#!/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.750579833984 57.130463513203 -gcp 10572.0 0.0 25.856151580811 57.141362263773 -gcp 10572.0 11715.0 25.878810882568 57.076574722762 -gcp 0.0 11715.0 25.774612426758 57.065190255237 LVVA_F7404_US1_GV898_0002.jpg LVVA_F7404_US1_GV898_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"507.0 8571.0\n522.0 7572.0\n513.0 354.0\n10062.0 447.0\n10023.0 10797.0\n2424.0 10674.0\n2490.0 9642.0\n2385.0 9657.0\n2220.0 9471.0\n1971.0 9351.0\n1903.5 9168.0\n1780.5 9033.0\n1615.5 9030.0\n1401.0 8991.0\n1216.5 8859.0\n1135.5 8763.0\n1036.5 8703.0\n874.5 8715.0\n714.0 8700.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV898_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV898_0002.vrt LVVA_F7404_US1_GV898_0002.cut.vrt

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