#!/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 316.7934278021638 756.5857496848707 20.96668067397325 56.49995701869416 -gcp 16054.888937716183 616.5960702822722 21.500071108394103 56.50000623516196 -gcp 16182.874599064146 13777.367359342887 21.499988857157597 56.24997850114886 -gcp 362.95054842520045 13916.16644747924 20.966636349695804 56.24997889874476 Lat075_006_Barta_ca1930.jpg Lat075_006_Barta_ca1930.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e \"315.72367101303917 752.2429596907785\n336.0 6863.170071783545\n64.49849548645935 6866.2682219768085\n72.00037612838514 9216.753313086692\n343.53159478435305 9214.512424075096\n362.0845035105316 13917.802595251243\n16184.813440320964 13779.028989508559\n16054.6165997994 615.0347874102706\" | \
  gdaltransform -tps -output_xy Lat075_006_Barta_ca1930.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}'))))" \
  Lat075_006_Barta_ca1930.vrt Lat075_006_Barta_ca1930.cut.vrt

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