#!/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 3858.0 3108.0 25.826666666666668 56 -gcp 23622.0 3120.0 26.276666666666667 56 -gcp 23721.0 23613.0 26.276666666666667 55.75 -gcp 3843.0 23550.0 25.826666666666668 55.75 lvva/Kurzemes,Vidzemes,Kaunas_un_Vitebskas_gubernas_Sederes_apkartne__6828_A1_L22_50_02.JPG Kurzemes,Vidzemes,Kaunas_un_Vitebskas_gubernas_Sederes_apkartne__6828_A1_L22_50_02.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "3858.0 3108.0\n23622.0 3120.0\n23718.0 18221.991671293727\n23540.25 18368.129650194336\n23283.0 18932.24125485841\n22926.0 19429.494169905607\n22416.0 20586.35313714603\n22269.0 21516.213214880623\n22252.5 21718.52248750694\n21993.0 22429.407551360353\n21873.0 22865.95252637424\n21829.5 23069.94503053859\n21654.0 23413.46751804553\n21583.5 23609.97001665741\n3843.0 23550.0" | \
  gdaltransform -tps -output_xy Kurzemes,Vidzemes,Kaunas_un_Vitebskas_gubernas_Sederes_apkartne__6828_A1_L22_50_02.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  Kurzemes,Vidzemes,Kaunas_un_Vitebskas_gubernas_Sederes_apkartne__6828_A1_L22_50_02.vrt Kurzemes,Vidzemes,Kaunas_un_Vitebskas_gubernas_Sederes_apkartne__6828_A1_L22_50_02.cut.vrt

# Generate tiles
gdal2tiles.py -r $GDAL_SAMPLING_TILE --xyz -z $GDAL_TILE_ZOOMS -x --processes=$GDAL_TILE_PROCESSES Kurzemes,Vidzemes,Kaunas_un_Vitebskas_gubernas_Sederes_apkartne__6828_A1_L22_50_02.cut.vrt Kurzemes,Vidzemes,Kaunas_un_Vitebskas_gubernas_Sederes_apkartne__6828_A1_L22_50_02.xyz