#!/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.273803710938 57.545036533239 -gcp 13860.0 0.0 26.420917510986 57.554523633391 -gcp 13860.0 12204.0 26.43928527832 57.490737584631 -gcp 0.0 12204.0 26.293201446533 57.479849636419 LVVA_F7404_US1_GV2080_0002.jpg LVVA_F7404_US1_GV2080_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"4101.0 525.0\n13266.0 465.0\n13284.0 4104.0\n12225.0 4161.0\n11769.0 4245.0\n11130.75 4394.25\n10996.5 4476.0\n10959.0 4415.25\n10850.25 4474.5\n10728.0 4439.25\n10410.75 4553.25\n10182.0 4414.5\n10143.0 4293.0\n9691.5 3927.0\n9495.0 3879.0\n9403.5 3967.5\n9568.5 4105.5\n9595.5 4191.0\n9504.0 4320.0\n9351.0 4402.5\n9154.5 4444.5\n9069.0 4479.0\n9036.0 4516.5\n9046.5 4578.0\n9099.0 4609.5\n9223.5 4563.0\n9282.0 4561.5\n9429.0 4699.5\n9516.0 4806.0\n9510.0 4862.25\n9470.25 4928.25\n9465.75 4962.0\n9618.0 5049.75\n9668.25 5085.0\n9706.5 5143.5\n9723.0 5226.0\n9703.5 5324.25\n9652.5 5405.25\n9510.75 5496.75\n9369.0 5568.75\n9298.5 5629.5\n9180.0 5515.5\n9133.5 5548.5\n9072.0 5505.0\n8916.0 5611.5\n8788.5 5458.5\n8689.5 5521.5\n8652.0 5496.0\n8646.0 5422.5\n8550.0 5427.0\n8500.5 5634.0\n8502.0 5742.0\n8413.5 5884.5\n8308.5 5895.0\n8281.5 5961.0\n8124.0 6007.5\n6783.0 6031.5\n4140.0 6046.5\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV2080_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV2080_0002.vrt LVVA_F7404_US1_GV2080_0002.cut.vrt

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