#!/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 24.813394546509 57.935540328129 -gcp 15690.0 0.0 24.829919469472 57.855466680092 -gcp 15690.0 20754.0 24.626369476318 57.840045604744 -gcp 0.0 20754.0 24.608087539673 57.919725722244 LVVA_F7404_US1_GV1938_0003.jpg LVVA_F7404_US1_GV1938_0003.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"3519.0 12240.0\n3513.0 2443.5\n3480.0 2022.0\n4863.0 696.0\n10594.5 7683.0\n10659.0 10881.0\n14886.75 11080.5\n14872.5 12306.0\n15315.0 12507.0\n15256.5 12774.0\n15207.0 12886.5\n15135.0 12997.5\n14853.0 13380.0\n14844.0 14661.0\n12243.0 14653.5\n12205.5 14910.0\n12096.0 15066.0\n11934.0 15115.5\n11838.0 15081.0\n11748.0 15030.0\n11556.0 14887.5\n11368.5 14778.0\n11272.5 14715.0\n11151.0 14709.0\n11011.5 14740.5\n10887.0 14800.5\n10774.5 14853.0\n10665.0 14956.5\n10555.5 15181.5\n10534.5 16009.5\n10603.5 16099.5\n10521.0 16551.0\n10462.5 19861.5\n3370.5 19891.5\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1938_0003.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1938_0003.vrt LVVA_F7404_US1_GV1938_0003.cut.vrt

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