#!/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.700210571289 57.368153468476 -gcp 18261.0 0.0 26.874103546143 57.38541305969 -gcp 18261.0 8223.0 26.886119842529 57.349635496397 -gcp 0.0 8223.0 26.709222793579 57.332451723974 LVVA_F7404_US1_GV1805_0002.jpg LVVA_F7404_US1_GV1805_0002.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"5710.5 2817.0\n5667.0 2787.0\n5560.5 2773.5\n5502.0 2697.0\n5409.0 2592.0\n5328.0 2506.5\n5197.5 2502.0\n5092.5 2559.0\n4939.5 2515.5\n4741.5 2322.0\n4735.5 2097.0\n4965.0 2148.0\n5140.5 2220.0\n5544.0 2383.5\n5586.0 2358.0\n5617.5 2308.5\n5694.0 2296.5\n5902.5 2482.5\n6052.5 2548.5\n6190.5 2629.5\n6321.0 2679.0\n6319.5 1131.0\n8877.0 420.0\n11211.0 429.0\n13101.0 435.0\n13062.0 5322.0\n9520.5 5310.0\n9492.0 7470.0\n5295.0 7500.0\n5304.0 6456.0\n6042.0 4878.0\n5844.0 4578.0\n5718.0 4470.0\n5629.5 4192.5\n5721.0 3952.5\n5728.5 3411.0\n5887.5 3375.0\n5964.0 3384.0\n6021.0 3428.25\n6077.25 3431.25\n6107.25 3399.0\n6113.25 3298.5\n6136.5 3267.0\n6177.0 3252.75\n6314.25 3245.25\n6306.0 2981.25\n6260.25 2944.5\n6241.5 2883.0\n6194.25 2916.0\n6120.75 2919.75\n6100.5 2907.75\n6087.75 2835.75\n6117.0 2813.25\n6069.0 2787.75\n6020.25 2771.25\n5966.25 2790.0\n5892.0 2830.5\n5776.5 2836.5\" | \
  gdaltransform -tps -output_xy LVVA_F7404_US1_GV1805_0002.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  LVVA_F7404_US1_GV1805_0002.vrt LVVA_F7404_US1_GV1805_0002.cut.vrt

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