#!/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 1227.0933279742767 567.8314131038312 35.99767807 64.00007586 -gcp 15238.878135048231 498.3364797334814 47.99795348 64.000257 -gcp 16261.57234726688 11066.344322598557 47.99820539 60.00018488 -gcp 284.6012875132635 11154.957304015168 35.99796387 60.00001024 -gcp 8233.875 856.9166435313715 41.99780374 64.00016927 -gcp 8278.875 11472.07745696835 41.99807408 60.00010031 -gcp 2389.5 652.7285067873304 36.99820635 63.9635592 -gcp 3556.5 727.7662141779789 37.99989987 63.93103544 -gcp 4723.5 789.3431372549019 38.99909497 63.90395074 -gcp 5898.0 828.3340874811463 40.0025178 63.88594487 -gcp 7066.5 846.3431372549019 40.99956077 63.87638718 -gcp 9403.5 843.341628959276 42.99279669 63.87294999 -gcp 10575.0 802.8212669683257 43.99315433 63.88690649 -gcp 11742.0 757.815233785822 44.990911 63.90268369 -gcp 12912.0 688.7805429864254 45.99374469 63.92813528 -gcp 14079.0 598.7450980392157 46.9972826 63.96207713 -gcp 15495.0 3141.101809954751 47.99833279 62.95046679 -gcp 15750.0 5783.024886877828 47.99785435 61.93449027 -gcp 16006.5 8424.0 47.997798 60.95318548 -gcp 14940.0 11187.461538461539 46.99781726 59.95992745 -gcp 13614.0 11286.511312217195 45.99798786 59.92753044 -gcp 12282.0 11365.91101055807 44.99667932 59.90211571 -gcp 10944.0 11418.437405731522 43.99339935 59.88620827 -gcp 9612.0 11460.0 42.99712235 59.87507612 -gcp 6945.0 11466.401206636501 41.00099244 59.87710687 -gcp 5613.0 11443.889894419306 40.00418624 59.88763361 -gcp 4275.0 11403.4419306184 39.00177554 59.90449656 -gcp 2941.5 11338.909502262442 38.0007248 59.92985413 -gcp 1609.91288833028 11256.800989923031 36.9986163 59.96143384 -gcp 519.0 8505.115384615385 35.99673969 60.95235946 -gcp 757.5 5862.0 36.00043156 61.93403598 -gcp 991.125 3214.9492835595775 35.99749594 62.95037096 -gcp 8256.0 6165.238310708899 41.99675399 61.80509342 -gcp 12011.25 6060.90346907994 44.9934242 61.83697725 -gcp 10668.0 3461.363499245852 43.9937073 62.83516264 -gcp 13089.0 3338.3016591251885 45.99614693 62.87828981 -gcp 13437.0 8638.024886877829 45.99582004 60.87846087 -gcp 10852.5 8766.582202111615 43.99386805 60.8366589 -gcp 4503.0 6097.93665158371 39.00356663 61.83783586 -gcp 5826.0 3486.0693815987934 40.0024449 62.83514129 -gcp 3400.5 3382.5173453996986 37.99836427 62.88025676 -gcp 3094.5 8685.23378582202 37.99991856 60.88153514 -gcp 5685.75 8795.489819004524 40.00508122 60.83626559 vlasenko/north/p-37_38.jpg p-37_38.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1227.0933279742767 567.8314131038312\n15238.878135048231 498.3364797334814\n16261.57234726688 11066.344322598557\n284.6012875132635 11154.957304015168" | \
  gdaltransform -tps -output_xy p-37_38.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  p-37_38.vrt p-37_38.cut.vrt

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