#!/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 2499.75 2551.6175653140635 25.333333333333332 56.25 -gcp 24279.75 2640.8450528071153 26.333333333333332 56.25 -gcp 24426.375 22760.00736520289 26.333333333333332 55.75 -gcp 2340.0 22669.81183991106 25.333333333333332 55.75 Latvijas_un_Lietuvas_topografiskas_kartes__Rokiski_un_apkartne__6828_A1_L32_7_02.JPG Latvijas_un_Lietuvas_topografiskas_kartes__Rokiski_un_apkartne__6828_A1_L32_7_02.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "13392.0 2661.036646307607\n2500.5 2549.2229317046085\n2433.0 12630.011660188782\n2340.0 22669.470849528043\n13381.5 22750.51998889506\n24426.0 22759.494169905607\n24354.0 12730.517490283175\n24282.0 2638.495835646863" | \
  gdaltransform -tps -output_xy Latvijas_un_Lietuvas_topografiskas_kartes__Rokiski_un_apkartne__6828_A1_L32_7_02.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  Latvijas_un_Lietuvas_topografiskas_kartes__Rokiski_un_apkartne__6828_A1_L32_7_02.vrt Latvijas_un_Lietuvas_topografiskas_kartes__Rokiski_un_apkartne__6828_A1_L32_7_02.cut.vrt

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