#!/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 27.013578414917 57.564606724948 -gcp 9126.0 0.0 27.028341293335 57.517159648319 -gcp 9126.0 19449.0 26.838226318359 57.500976675073 -gcp 0.0 19449.0 26.825523376465 57.548260582254 LVVA_F7404_US1_GV1771_0002.jpg LVVA_F7404_US1_GV1771_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"856.5 9433.5\n858.0 4554.0\n7869.0 4614.0\n7861.5 6679.5\n7833.0 6679.5\n7693.5 6568.5\n7527.0 6460.5\n7318.5 6358.5\n7060.5 6312.0\n6870.0 6313.5\n6844.5 6408.0\n6807.0 6474.75\n6747.75 6505.5\n6544.5 6535.5\n6465.0 6975.0\n6479.25 7119.75\n6612.0 7098.0\n6666.75 7200.75\n6723.0 7268.25\n6763.5 7563.75\n6886.5 7870.5\n7016.25 8050.5\n7003.5 8220.0\n7012.5 8389.5\n6642.75 9501.0\n6686.25 9503.25\n6808.5 9572.25\n6855.0 9618.75\n6924.0 9687.0\n6978.75 9651.0\n6995.25 9573.75\n7058.25 9555.75\n7192.5 9594.0\n7335.0 9598.5\n7687.5 9528.0\n7746.0 9475.5\n7792.5 9345.0\n7836.0 9301.5\n7833.75 9704.25\n7869.0 9709.5\n7920.0 9682.5\n8001.0 9687.75\n8073.75 9721.5\n8083.5 9904.5\n8160.0 9939.0\n8214.0 10171.5\n8250.0 10432.5\n8188.5 10561.5\n8008.5 10728.75\n7842.0 10731.75\n7852.5 12540.0\n7938.0 12598.5\n7927.5 12957.0\n7968.0 13152.0\n7984.5 13738.5\n7992.0 18522.0\n894.0 18495.0\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1771_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1771_0002.vrt LVVA_F7404_US1_GV1771_0002.cut.vrt

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