#!/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.037927627563 57.404606545616 -gcp 11676.0 0.0 25.064449310303 57.346092869791 -gcp 11676.0 20613.0 24.856910705566 57.328629441855 -gcp 0.0 20613.0 24.842147827148 57.386569578321 LVVA_F7404_US1_GV198_0002.jpg LVVA_F7404_US1_GV198_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"3954.0 9498.0\n613.5 9481.5\n708.0 5653.5\n4015.5 5827.5\n6324.0 5787.0\n8355.0 5760.0\n8682.0 6702.0\n8187.0 8407.5\n7560.0 9133.5\n7737.0 9327.0\n8026.5 9352.5\n8308.5 9633.0\n8469.0 9717.0\n8547.0 9813.0\n8449.5 9952.5\n8541.0 10101.0\n9048.0 10446.0\n9198.0 10623.0\n9352.5 10620.0\n9531.0 10827.0\n9418.5 11077.5\n9436.5 11374.5\n9384.0 11541.0\n9202.5 11700.0\n9085.5 11901.0\n8889.0 12087.0\n8848.5 12291.0\n8772.0 12475.5\n8746.5 12649.5\n8655.0 12843.0\n8619.0 13075.5\n8740.5 13266.0\n9046.5 13243.5\n9328.5 13083.0\n9564.0 13198.5\n9730.5 12759.0\n9894.0 12768.0\n10014.0 12844.5\n9924.0 12996.0\n9936.0 13107.0\n10132.5 13327.5\n10195.5 13524.0\n10069.5 13605.0\n10131.0 13779.0\n10053.0 14047.5\n9900.0 14122.5\n9979.5 14413.5\n9858.0 14865.0\n9894.0 15066.0\n10044.0 15055.5\n10146.0 14958.0\n10284.0 15058.5\n10387.5 15045.0\n10366.5 14844.0\n10512.0 14835.0\n10593.0 14917.5\n10834.5 15015.0\n11055.0 14952.0\n11052.0 19788.0\n3900.0 19674.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV198_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV198_0002.vrt LVVA_F7404_US1_GV198_0002.cut.vrt

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