#!/usr/bin/env bash

geoXYZ="+proj=geocent +ellps=WGS84 +datum=WGS84 +no_defs"
isn93='+proj=lcc +lat_1=64.25 +lat_2=65.75 +lat_0=65 +lon_0=-19 +x_0=500000 +y_0=500000 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
isn2016='+proj=lcc +lat_1=64.25 +lat_2=65.75 +lat_0=65 +lon_0=-19 +x_0=2700000 +y_0=300000 +ellps=GRS80 +units=m +no_defs'
utm26n='+proj=utm +zone=26 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
utm27n='+proj=utm +zone=27 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
utm28n='+proj=utm +zone=28 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
lonlat='+proj=longlat +datum=WGS84 +no_defs'
latlon='+proj=latlong +datum=WGS84 +no_defs'

lzwparam='-co compress=lzw -co tiled=yes -co blockxsize=512 -co blockysize=512 -co predictor=2'
lzwparampy='--co compress=lzw --co tiled=yes --co blockxsize=512 --co blockysize=512 --co predictor=2'
demzstd='-co compress=zstd -co tiled=yes -co blockxsize=512 -co blockysize=512 -co predictor=2 -co zstd_level=15 -co num_threads=32'
imgzstd='-co compress=zstd -co tiled=yes -co bigtiff=yes -co blockxsize=1024 -co blockysize=1024 -co predictor=2 -co zstd_level=1 -co num_threads=32'
demzstdpy='--co compress=zstd --co tiled=yes --co blockxsize=512 --co blockysize=512 --co predictor=2 --co zstd_level=15 --co num_threads=32'

geoid2016=/u01/jmcb/1geodesy/geoid/Icegeoid_ISN2016.gtx
geoid93=/u01/jmcb/1geodesy/geoid/Icegeoid_ISN93.gtx

get_georef()
{
    xymin=$(gdalinfo $1 | grep 'Lower Left' | sed -e 's/[,()]//g' | awk '{print $3,$4}')
    xymax=$(gdalinfo $1 | grep 'Upper Right' | sed -e 's/[,()]//g' | awk '{print $3,$4}')
    echo $xymin $xymax
}

get_ullr()
{
    gdalinfo -nomd -norat $1 \
        | egrep -e 'Upper Left|Lower Right' \
        | sed -e 's/[,()]//g' \
        | sed -e 's/-/ -/g' \
        | awk '{print $3,$4,$10,$11}' \
        | tr '\n' ' '
}

get_px()
{
    gdalinfo -nomd -norat $1 \
        | grep 'Pixel Size' \
        | sed -e 's/)//g' \
        | awk -F ",-" '{print $2}'
}

get_ty()
{
    gdalinfo -nomd -norat $1 \
        | grep 'Type' \
        | awk '{print $4}' \
        | sed -e 's/[,()]//g' \
        | awk -F "=" '{print $2}'
}

get_proj()
{
    gdalsrsinfo --single-line -o proj4 $1
}

get_center()
{
    gdalinfo $1 \
    | grep 'Size is' \
    | sed 's/,//g' \
    | awk -v k=$2 '{print k*$3/2,k*$4/2}'
}

orto2jp2()
{
    nam=$(basename $1 .tif)
    gdal_translate -of JP2OpenJPEG -co INSPIRE_TG=YES -a_nodata 0 $nam.tif $nam.jp2
}

pl2zip1zip()
{
    unzip '*zip' 
    zip -r $1 6* 
    rm -r 6*
}

view_shade()
{
#https://github.com/dshean/imview
    mos_fn=$1
    gdaldem hillshade $mos_fn ${mos_fn}_shade.tif
    gdaladdo -ro -r average --config COMPRESS_OVERVIEW LZW --config BIGTIFF_OVERVIEW YES $mos_fn 2 4 8 16 32 64
    gdaladdo -ro -r average --config COMPRESS_OVERVIEW LZW --config BIGTIFF_OVERVIEW YES ${mos_fn}_shade.tif 2 4 8 16 32 64
    imviewer.py $mos_fn -overlay ${mos_fn}_shade.tif -label 'Elevation (m WGS84)' -of png -dpi 300 -scale x
}

view_fancyshade()
{
#https://github.com/dshean/imview
    mos_fn=$1
    gdal_translate -co compress=lzw $1 ~/cal/tmpnozip.tif

    whitebox_tools -r=TimeInDaylight -v --az_fraction=15.0 --max_dist=100.0 --lat=65 --long=-19 -i ~/cal/tmpnozip.tif -o ${mos_fn}_fancyshade.tif
    gdaladdo -ro -r average --config COMPRESS_OVERVIEW LZW --config BIGTIFF_OVERVIEW YES $mos_fn 2 4 8 16 32 64
    gdaladdo -ro -r average --config COMPRESS_OVERVIEW LZW --config BIGTIFF_OVERVIEW YES ${mos_fn}_fancyshade.tif 2 4 8 16 32 64
    imviewer.py $mos_fn -overlay ${mos_fn}_fancyshade.tif -label 'Elevation (m WGS84)'
}

geo2ell_isn2016()
{
    #nam="${1%.*}"
    nam=$(basename $1 .tif)
    echo $nam
    r=$(dirname "$1")
    gdalwarp -overwrite -r cubic -te $(get_georef $1) -tr $(get_px $1) $(get_px $1) \
        -t_srs "$(get_proj $1)" $lzwparam -of GTiff $geoid2016 ~/cal/geoid_frm_${nam}.tif
    gdal_calc.py --overwrite -A $1 -B ~/cal/geoid_frm_${nam}.tif \
        --outfile $dir/${nam}_zmae.tif --calc="A+B" $demzstdpy --type Float32
    rm ~/cal/geoid_frm_${nam}.tif
}

ell2geo_isn2016()
{
    #nam="${1%.*}"
    nam=$(basename $1 .tif)
    dir=$(dirname "$1")
    gdalwarp -overwrite -r cubic -te $(get_georef $1) -tr $(get_px $1) $(get_px $1) \
        -t_srs "$(get_proj $1)" $lzwparam -of GTiff $geoid2016 ~/cal/geoid_frm_${nam}.tif
    gdalinfo ~/cal/geoid_frm_${nam}.tif
    gdal_calc.py --overwrite -A $1 -B ~/cal/geoid_frm_${nam}.tif \
        --outfile $dir/${nam}_zmasl.tif --calc="A-B" $demzstdpy --type Float32
    rm ~/cal/geoid_frm_${nam}.tif
    #Faster tools, but not supporting zstd:
    #otbcli_BandMath -il $1 ~/cal/geoid_frm_${NAM}.tif -out $DIR/${NAM}_zmae.tif -exp "im1 - im2" 
    #image_calc -c "var_0 - var_1" $1 ~/cal/geoid_frm_${NAM}.tif -o $DIR/${NAM}_zmae.tif -d float32
}

geo2ell_isn93()
{
    nam=$(basename "$1" .tif)
    dir=$(dirname "$1")
    gdalwarp -overwrite -r cubic -te $(get_georef $1) -tr $(get_px $1) $(get_px $1) \
        -t_srs "$(get_proj $1)" $lzwparam -of GTiff $geoid93 ~/cal/geoid_frm_${nam}.tif
    gdal_calc.py --overwrite -A $1 -B ~/cal/geoid_frm_${nam}.tif \
        --outfile $dir/${nam}_zmae.tif --calc="A+B" $demzstdpy --type Float32
    rm ~/cal/geoid_frm_${nam}.tif
}

ell2geo_isn93()
{
    nam="${1%.*}"
    dir=$(dirname "$1")
    gdalwarp -overwrite -r cubic -te $(get_georef $1) -tr $(get_px $1) $(get_px $1) \
        -t_srs "$(get_proj $1)" $lzwparam -co BIGTIFF=YES -of GTiff $geoid93 ~/cal/geoid_frm_${nam}.tif
    gdal_calc.py --overwrite -A $1 -B ~/cal/geoid_frm_${nam}.tif \
        --outfile $dir/${nam}_zmasl.tif --calc="A-B" $demzstdpy --co BIGTIFF=YES --type Float32
    rm ~/cal/geoid_frm_${nam}.tif
    #Faster tools, but not supporting zstd:
    #otbcli_BandMath -il $1 ~/cal/geoid_frm_${NAM}.tif -out $DIR/${NAM}_zmae.tif -exp "im1 - im2" 
    #image_calc -c "var_0 - var_1" $1 ~/cal/geoid_frm_${NAM}.tif -o $DIR/${NAM}_zmae.tif -d float32
}

shiftXYZ()
{
    nam="${1%.*}"
    dir=$(dirname "$1")
    
    #x_rou=$(echo $2 | printf "%.3f" )
    #y_rou=$(printf "%.3f" $3)
    #z_rou=$(printf "%.3f" $4)
    
    ULLR=$(get_georef $1 | awk -v dX=$2 -v dY=$3 '{ printf "%.3f %.3f %.3f %.3f", $1+dX, $4+dY, $3+dX, $2+dY}')
    gdal_translate -of VRT -a_ullr $ULLR $1 -a_srs "$(get_proj $1)" ~/cal/${nam}_sftXY.vrt
    image_calc -c "var_0+$4" --tile-size 512 512 -d float32 ~/cal/${nam}_sftXY.vrt -o $dir/${nam}_x${2}_y${3}_z${4}.tif
    rm ~/cal/${nam}_sftXY.vrt
}

shiftXY()
{
    nam="${1%.*}"
    dir=$(dirname "$1")
    
    x_rou=$(printf "%.2f" $2)
    y_rou=$(printf "%.2f" $3)
    
    ULLR=$(get_georef $1 | awk -v dX=$2 -v dY=$3 '{ printf "%.3f %.3f %.3f %.3f", $1+dX, $4+dY, $3+dX, $2+dY}')
    gdal_translate $imgzstd -a_ullr $ULLR -a_srs "$(get_proj $1)" $1 ${nam}_sft_x${x_rou}_y${y_rou}.tif
}

################################################################
#################### DEBUG & OLD STUFF ZONE #################### 
################################################################

#geo2ell()
#{
#    NAM=$(basename "$1" | sed 's/.\{4\}$//')
#    DIR=$(dirname "$1")
#    gdal_translate -of VRT -a_srs "${isn93}" $1 ~/cal/${NAM}.vrt
#    gdalwarp -overwrite -r bilinear -te $(get_georef $1) -tr $(get_px $1) $(get_px $1) -t_srs "${isn93}" -of GTiff $geoid ~/cal/geoid-frm-${NAM}.tif
#    image_calc -c "var_0 + var_1" --tile-size 512,512 ~/cal/${NAM}.vrt ~/cal/geoid-frm-${NAM}.tif -o $DIR/${NAM}-zmae.tif -d float32 $2
#    rm ~/cal/${NAM}.vrt ~/cal/geoid-frm-${NAM}.tif
#}
#
#ell2geo()
#{
#    NAM=$(basename "$1" | sed 's/.\{4\}$//')
#    DIR=$(dirname "$1")
#    gdal_translate -of VRT -a_srs "${isn2016}" $1 ~/cal/${NAM}.vrt
#    gdalwarp -overwrite -r bilinear -te $(get_georef $1) -tr $(get_px $1) $(get_px $1) -t_srs "${isn2016}" -of GTiff $geoid ~/cal/geoid-frm-${NAM}.tif
#    gdal_calc.py --overwrite -A ~/cal/${NAM}.vrt -B ~/cal/geoid-frm-${NAM}.tif --outfile $DIR/${NAM}-zmasl.tif --calc="A-B" $demzstdpy --type Float32
#    #image_calc -c "var_0 - var_1" --tile-size 512,512 ~/cal/${NAM}.vrt ~/cal/geoid-frm-${NAM}.tif -o $DIR/${NAM}-zmasl.tif -d float32 $2
#    rm ~/cal/${NAM}.vrt ~/cal/geoid-frm-${NAM}.tif
#}