#!/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 2895.2722014173887 4477.8798329654155 24.05246124 56.94126239 -gcp 21498.0000022817 3.6445378896426135e-06 24.07461512 56.94363214 -gcp 20619.138114408648 19206.47968930214 24.07299489 56.93163949 -gcp 2976.7884915758295 19869.302913988915 24.05231774 56.93132758 GX4799/GX_4799_SD_frame_111.jpg GX_4799_SD_frame_111.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "2616.0 20154.003342618387\n2190.0 19668.066852367687\n930.0 19662.157103064066\n852.0 2447.9665738161566\n2628.0 2484.2473537604455\n2652.0 1128.441225626741\n18012.0 960.3275766016714\n18000.0 3025.490807799443\n21288.0 2989.4373259052927\n21216.0 21359.959888579386\n2664.0 21431.993314763233" | \
  gdaltransform -tps -output_xy GX_4799_SD_frame_111.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  GX_4799_SD_frame_111.vrt GX_4799_SD_frame_111.cut.vrt

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