#!/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 4006.5 1231.6734296831573 23.1266666 57.25 -gcp 33534.825694264706 1131.72410191569 24.4766666 57.25 -gcp 33829.5 30693.167176209005 24.4766666 56.5 -gcp 3706.1759500078365 30787.665370974268 23.1266666 56.5 Latvijas_un_Igaunijas_topografiskas_kartes_Jelgava_6828_A1_L36_19_02.JPG Latvijas_un_Igaunijas_topografiskas_kartes_Jelgava_6828_A1_L36_19_02.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "18768.0 1213.5183231538035\n4005.0 1230.0\n3804.0 16011.001665741254\n3702.75 30785.247917823435\n18681.0 30840.056635202665\n33831.0 30693.0\n33687.0 15913.500832870628\n33534.0 1131.0" | \
  gdaltransform -tps -output_xy Latvijas_un_Igaunijas_topografiskas_kartes_Jelgava_6828_A1_L36_19_02.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  Latvijas_un_Igaunijas_topografiskas_kartes_Jelgava_6828_A1_L36_19_02.vrt Latvijas_un_Igaunijas_topografiskas_kartes_Jelgava_6828_A1_L36_19_02.cut.vrt

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