#!/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 1310.25 350.4952803997779 47.99795348 64.000257 -gcp 15045.0 462.9998611882288 59.99831835 64.00041159 -gcp 15918.75 10893.202803997778 59.99852535 60.00033392 -gcp 255.75 10761.952387562465 47.99820539 60.00018488 -gcp 8177.0625 738.2348695169351 53.99812565 64.00033807 -gcp 8090.0625 11189.284355913382 53.99835636 60.00026304 -gcp 2451.0 459.27903469079934 49.00232539 63.95958107 -gcp 3594.0 552.3257918552036 50.00496875 63.92539204 -gcp 4740.0 630.0 51.00975326 63.89799701 -gcp 5887.5 679.8898944193062 52.00789871 63.88101666 -gcp 7038.0 718.9095022624435 53.00985673 63.86896131 -gcp 9327.0 739.765460030166 55.00117064 63.86813354 -gcp 10473.0 717.2541478129713 55.99866432 63.88115813 -gcp 11616.0 684.2699849170436 56.99470136 63.89845105 -gcp 12765.0 628.7420814479638 57.99810082 63.92497867 -gcp 13909.5 555.2051282051282 59.00036008 63.95887326 -gcp 15267.0 3069.9072398190046 60.0014293 62.95103527 -gcp 15484.5 5676.94532428356 60.00065176 61.93515805 -gcp 15700.5 8281.673453996984 59.99889033 60.95298059 -gcp 14623.5 10992.871794871797 59.00062014 59.96067117 -gcp 13321.5 11073.912518853696 58.00086304 59.92771198 -gcp 12013.5 11131.944947209653 56.9993022 59.90300262 -gcp 10705.5 11172.428355957767 55.99961769 59.88458789 -gcp 9399.0 11191.919306184012 55.00234929 59.8736911 -gcp 6784.5 11166.381598793363 53.00749216 59.87495762 -gcp 5476.5 11122.865007541479 52.00843422 59.88658231 -gcp 4168.5 11062.85294117647 51.0078736 59.90411459 -gcp 2863.5 10981.90346907994 50.0072775 59.92916594 -gcp 1555.5 10875.349924585218 49.00100824 59.96341246 -gcp 520.5 8152.7209653092 47.99852761 60.95374519 -gcp 780.0 5553.439668174962 47.99540963 61.93564753 -gcp 1044.0 2952.4208144796376 47.99572815 62.95100563 -gcp 8142.0 5959.7375565610855 54.00923838 61.80354038 -gcp 11823.783121058586 5907.279725422555 57.00452985 61.83529117 -gcp 12912.0 3240.0542986425335 58.00549344 62.87597065 -gcp 10537.5 3330.099547511312 56.00430534 62.83262931 -gcp 13189.477431283625 8460.401186719047 58.00593342 60.87894785 -gcp 10652.977431283625 8559.450960474702 56.00379989 60.83445692 -gcp 4461.75 5843.491704374057 51.0135274 61.83626906 -gcp 5793.0 3294.330316742082 52.01510355 62.83125805 -gcp 3414.0 3160.758672699849 50.00787549 62.87626998 -gcp 5587.5 8513.086726998492 52.01483832 60.83593267 -gcp 3051.0 8373.343137254902 50.01053363 60.8797891 vlasenko/north/p-39_40.jpg p-39_40.vrt

# Add cutline to VRT
gdalwarp -r $GDAL_SAMPLING_WARP -tps -dstalpha -cutline_srs EPSG:4326 \
  -cutline "POLYGON(($(echo -e "1310.25 350.4952803997779\n15045.0 462.9998611882288\n15918.75 10893.202803997778\n255.75 10761.952387562465" | \
  gdaltransform -tps -output_xy p-39_40.vrt | \
  awk 'NR==1{first=$0} {printf "%s %s,", $1,$2} END{print " " first}')))" \
  p-39_40.vrt p-39_40.cut.vrt

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