#!/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.017242431641 57.574963158321 -gcp 15060.0 0.0 25.170879364014 57.587064943843 -gcp 15060.0 16833.0 25.192337036133 57.500238818634 -gcp 0.0 16833.0 25.047969818115 57.486954851951 LVVA_F7404_US1_GV1911_0002.jpg LVVA_F7404_US1_GV1911_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"765.0 8562.0\n318.0 8298.0\n306.0 4038.0\n741.0 3714.0\n735.0 1536.0\n2670.0 1548.0\n5724.0 3222.0\n6156.0 2760.0\n7068.0 3135.0\n7914.0 3552.0\n9315.0 3675.0\n9510.0 3492.0\n9714.0 3429.0\n14235.0 3426.0\n14439.0 15615.0\n822.0 15642.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1911_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1911_0002.vrt LVVA_F7404_US1_GV1911_0002.cut.vrt

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