#!/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 26.01413268195 57.016152606719 -gcp 10479.0 0.0 26.112770288182 57.026282998371 -gcp 10479.0 8676.0 26.126358284315 56.982168346255 -gcp 0.0 8676.0 26.027720678083 56.972037954603 LVVA_F7404_US1_GV949_0002.jpg LVVA_F7404_US1_GV949_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"558.0 3978.0\n1572.0 3759.0\n2124.0 3390.0\n2601.0 3393.0\n2796.0 3438.0\n3063.0 3609.0\n3237.0 3630.0\n3378.0 3720.0\n3444.0 3696.0\n3663.0 3198.0\n3784.5 3058.5\n3967.5 2979.0\n4170.0 2925.0\n4363.5 2968.5\n4486.5 3024.0\n4528.5 3088.5\n4540.5 3198.0\n4522.5 3312.0\n4459.5 3445.5\n4416.0 3513.0\n4264.5 3706.5\n4225.5 3775.5\n4281.0 3808.5\n4375.5 3955.5\n4446.0 3943.5\n4669.5 3769.5\n4801.5 3595.5\n4860.0 3429.0\n4884.0 3277.5\n5712.0 3210.0\n6765.0 2691.0\n6486.0 2031.0\n6666.0 1770.0\n7722.0 1443.0\n10020.0 2097.0\n10011.0 7671.0\n537.0 7674.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV949_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV949_0002.vrt LVVA_F7404_US1_GV949_0002.cut.vrt

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