#!/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.588445663452 57.228390340745 -gcp 14847.0 0.0 25.729637145996 57.241628902966 -gcp 14847.0 8898.0 25.742168426514 57.196086898356 -gcp 0.0 8898.0 25.603551864624 57.18301805405 LVVA_F7404_US1_GV1572_0002.jpg LVVA_F7404_US1_GV1572_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"714.0 549.0\n10062.0 555.0\n10083.0 4239.0\n8328.0 5935.5\n6294.0 5610.0\n6266.25 5462.25\n6041.25 5539.5\n5820.0 5445.75\n5625.0 5440.5\n5596.5 5566.5\n5523.0 5682.0\n5461.5 5808.0\n5500.5 5880.0\n5503.5 6009.0\n5614.5 6135.0\n5739.0 6243.0\n5628.0 6349.5\n5604.0 6555.0\n5487.0 6631.5\n5364.0 6819.0\n5206.5 7039.5\n4929.0 7303.5\n4692.0 7383.0\n4567.5 7428.0\n4440.0 7392.0\n4345.5 7306.5\n4282.5 7161.0\n4096.5 7021.5\n4030.5 6921.0\n3982.5 6793.5\n3967.5 6685.5\n3906.0 6628.5\n3817.5 6621.0\n3169.5 6460.5\n2436.0 6747.0\n1893.0 6822.0\n687.0 6852.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1572_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1572_0002.vrt LVVA_F7404_US1_GV1572_0002.cut.vrt

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