#!/bin/sh
# # Animation : sea-ice cover in Arctic sea (2016).
# ## Stereographic projection (-Js) with some regions leaking out of the frame
#   based on (IBCAO map)[https://www.ngdc.noaa.gov/mgg/bathymetry/arctic/AGU_Fall_2011_IBCAO.pdf]
#
# ## Dataset :
#   (topography and bathymetry)[https://docs.generic-mapping-tools.org/dev/datasets/remote-data.html?highlight=remote%20data#global-earth-relief-grids]
#   (sea-ice concentration)[https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-sea-ice-concentration?tab=overview]
#
#   The dataset are meant to be read from ${work_dir} such as :
#       ${work_dir}/north/[...] = north hemisphere
#       ${work_dir}/south/[...] = south hemisphere 
#   (see include and pre.sh)
#
# ## Layers :
#    The goal is to have the ocean and the lands confined in an inner area with some extra lands leaking outside
#
# (1) INNER             = Ocean & lands (@earth_relief)                                 => PRE.sh
# (2) INNER/[OUTTER]    = Ice (dataset) / [shifted (-dx -dy) and trimmed (clip)]        => MAIN.sh
# (3) INNER             = Basemap / [re-shifted (+dx +dy)]                              => POST.sh
# (4) OUTTER            = Land (@earth_relief) shifted (-dx -dy) and trimmed (coast -E) => POST.sh
# (5) colorbar      
#
# ## Functions :
#   To substentially improve code readability (fewer repeated lines), I resorted to use function
#
# (1) gmt_project   = defines the projection, scales and calculates the shift between INNER and OUTTER regions
# (2) gmt_plotb     = background layers (topography)
# (3) gmt_plotm     = main layers (sea-ice) + title
# (4) gmt_plotf     = foreground layers (extra pieces of land) + colorbars
#
# (nota) : subtitle calls LaTeX routines
#
# ### Credits :
# Guillaume D.-A. 2021-09 (version 2)
# # # # # # # # # # # # # # # # # # #


# # # # # # # # # # # # # # # # # # # # # # # # # #
#             DOWNLOAD API (python)               #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# /!\ Requires an account on (ECMWF copernicus)[https://cds.climate.copernicus.eu]

cat << 'EOF' > download.py  
import cdsapi
import argparse


    parser = argparse.ArgumentParser()
    parser.add_argument("number", type=int)
    args = parser.parse_args()

    c = cdsapi.Client()

    c.retrieve(
    'satellite-sea-ice-concentration',
    {
        'origin': 'esa_cci',
        'region': 'northern_hemisphere','southern_hemisphere',
        'cdr_type': 'cdr',
        'year': str(args.number),
        'month': [
            '01', '02', '03',
            '04', '05', '06',
            '07', '08', '09',
            '10', '11', '12',
        ],
        'day': [
            '01', '02', '03',
            '04', '05', '06',
            '07', '08', '09',
            '10', '11', '12',
            '13', '14', '15',
            '16', '17', '18',
            '19', '20', '21',
            '22', '23', '24',
            '25', '26', '27',
            '28', '29', '30',
            '31',
        ],
        'version': 'v2',
        'variable': 'all',
        'format': 'zip',
    },
    "download_" + str(args.number) + ".zip")
EOF

# # # # # # # # # # # # # # # # # # # # # # # # # #
#             DOWNLOAD API (python)               #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
#                 FUNCTIONS\                      #
# # # # # # # # # # # # # # # # # # # # # # # # # #

cat << 'EOF' > gmt_project.sh
#!/bin/sh
#    # longitude center
#    # latitude center
#    # latitude inner circle
#    # latitude outter circle
#    # latitude true scale
#    # scale

    lon_ref=$1
    lat_ref=$2

    lat_in=$3
    lat_out=$4
    lat_scl=$5
    scale="1:${6}"

    PROJECTION="-Js${lon_ref}/${lat_ref}/${lat_scl}/${scale}"

    # Sort latitude order depending on North or South hemisphere
    if [ $lat_ref -gt $lat_in ];
    then
        INNER="-R0/360/${lat_in}/${lat_ref}"
        OUTTER="-R0/360/${lat_out}/${lat_ref}"
    else
        INNER="-R0/360/${lat_ref}/${lat_in}"
        OUTTER="-R0/360/${lat_ref}/${lat_out}"
    fi

    # Get the dimensions of map (inner and outter) 
    W1=$(gmt mapproject ${INNER} ${PROJECTION}  -Ww)
    H1=$(gmt mapproject ${INNER} ${PROJECTION} -Wh)
    W2=$(gmt mapproject ${OUTTER} ${PROJECTION} -Ww)
    H2=$(gmt mapproject ${OUTTER} ${PROJECTION} -Wh)

    # Calculate the shift required to match inner and outter
    dx=$(gmt math -Q ${W1} ${W2} SUB 2 DIV =)
    dy=$(gmt math -Q ${H1} ${H2} SUB 2 DIV =)

    # output to be read by the script
    echo ${PROJECTION}
    echo ${INNER}
    echo ${OUTTER}
    echo ${dx} 
    echo ${dy}
    echo ${H1}  # used for colorbar size [ice]
    echo ${W1}  # used for colorbar size [topo]
EOF

cat << 'EOF' > gmt_plotb.sh
#!/bin/sh
#    # Projection
#    # Region
#    # Map shift x
#    # Map shift y

    projection=$1
    region=$2
    shift_x=$3
    shift_y=$4


    gmt makecpt -Cbukavu -T-4500/1500 -H > relief.cpt
    gmt makecpt -Cvik -G-0.5/0 -T0/101/1 -H -N > ice.cpt

    gmt grdimage ${projection} ${region} @earth_relief_01m -Crelief.cpt -I+d -X${shift_x}c -Y${shift_y}c
EOF

cat << 'EOF' > gmt_plotm.sh
#!/bin/sh
#    # Projection
#    # Region
#    # Map shift x
#    # Map shift y
#    # Data

    projection=$1
    region=$2
    shift_x=$3
    shift_y=$4
    data=$5

    # Reproject the data from distance [in km] from the pole to [lon / lat] and oversample using bilinear interpolation
    gmt grdproject ${projection} ${region} ${data} -Gdata_bis.nc -C -I -Fk
    gmt grdsample data_bis.nc+n0 -I1m -nl -Gdata_ter.nc

    gmt grdimage ${projection} ${region} data_ter.nc -Cice.cpt -Q -X${shift_x} -Y${shift_y}
EOF

cat << 'EOF' > gmt_plotf.sh
#!/bin/sh
#    # Projection
#    # Region
#    # Map shift x
#    # Map shift y
    
    projection=$1
    region=$2
    shift_x=$3
    shift_y=$4
    
    gmt basemap ${projection} ${region} -Ba60f10 -X${shift_x} -Y${shift_y}
    
EOF

# # # # # # # # # # # # # # # # # # # # # # # # # #
#                 /FUNCTIONS                      #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
#               INCLUDE.SH\                       #
# # # # # # # # # # # # # # # # # # # # # # # # # #

cat << 'EOF' > include.sh
    work_dir="/home/user/gmt_poles"
    
    DIR_north="${work_dir}/north"
    DIR_south="${work_dir}/south"

    # Pull the functions made "on the fly" in working directory
    cp ${work_dir}/gmt_p*.sh ./
    chmod +x ./gmt_p*.sh

    lat_in=65       # Inner circle extent
    lat_out=50      # Outter circule extent
    lat_scl=77      # Latitude of true scaling
    scale=93300000  # 1:scale <=> here roughly 6cm wide
    shift=12        # offset between South pole map and North pole map
    x_origin=3      # offset of the North pole map to have a quasi centered picture
    y_origin=3      # offset of the North pole map to have the colorbar visible [topo]

    # # # # NORTH HEMISPHERE
    {
    read -r PROJECTION_north
    read -r INNER_north
    read -r OUTTER_north
    read -r dx_north
    read -r dy_north
    read -r HEIGHT_map
    read -r WIDTH_map
    } <<< "$(./gmt_project.sh 0 90 ${lat_in} ${lat_out} ${lat_scl} ${scale})"

    # # # # # SOUTH HEMISPHERE
    {
    read -r PROJECTION_south
    read -r INNER_south
    read -r OUTTER_south
    read -r dx_south
    read -r dy_south
    } <<< "$(./gmt_project.sh 0 -90 -${lat_in} -${lat_out} -${lat_scl} ${scale})"

EOF

# # # # # # # # # # # # # # # # # # # # # # # # # #
#               /INCLUDE.SH                       #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
#                   PRE.SH\                       #
# # # # # # # # # # # # # # # # # # # # # # # # # #

cat << 'EOF' > pre.sh
gmt begin

    # List of files to process
    ls -v ${DIR_north}/north_* | awk -F[_.] '{print $2}' > list.txt


    # # # # NORTH HEMISPHERE
    ./gmt_plotb.sh ${PROJECTION_north} ${INNER_north} ${x_origin} ${y_origin}


    # # # # # SOUTH HEMISPHERE
    ./gmt_plotb.sh ${PROJECTION_south} ${INNER_south} ${shift} 0

gmt end 
EOF

# # # # # # # # # # # # # # # # # # # # # # # # # #
#                   /PRE.SH                       #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
#                  MAIN.SH\                       #
# # # # # # # # # # # # # # # # # # # # # # # # # #

cat << 'EOF' > main.sh
gmt begin

    # # # # NORTH HEMISPHERE
    ../gmt_plotm.sh ${PROJECTION_north} ${INNER_north} ${x_origin} ${y_origin} ${DIR_north}/north_${MOVIE_COL0}.nc?"ice_conc"
    

    # # # # SOUTH HEMISPHERE
    ../gmt_plotm.sh ${PROJECTION_south} ${OUTTER_south} $(gmt math -Q ${shift} ${dx_south} ADD =) ${dy_south} ${DIR_south}/south_${MOVIE_COL0}.nc?"ice_conc"
    
    # Title (+subtitle using LaTeX) positioning
    shf_t=$(gmt math -Q -fT ${MOVIE_COL0}T --FORMAT_DATE_IN=yyyymmdd --FORMAT_DATE_OUT=yyyy-o-dd = | awk -F[-] '{print $2,$1}')
    shf_x=$(gmt math -Q ${dx_south} -1 MUL ${shift} 2 DIV SUB =)
    shf_y=$(gmt math -Q ${dy_south} -1.5 MUL =)

    gmt basemap ${PROJECTION_south} ${INNER_south} -B+t"${shf_t}"+s"@[\textsc{sea-ice concentration}@[" -X${shf_x} -Y${shf_y}

gmt end
EOF

# # # # # # # # # # # # # # # # # # # # # # # # # #
#                  /MAIN.SH                       #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
#                  POST.SH\                       #
# # # # # # # # # # # # # # # # # # # # # # # # # #

cat << 'EOF' > post.sh
gmt begin

    # # # # NORTH HEMISPHERE
    ./gmt_plotf.sh ${PROJECTION_north} ${INNER_north} $x_origin $y_origin

    # Add the Nunavut islands (only, clip out the mainland province), Greenland and Iceland
    gmt coast -ECA.NU,GL,IS+c ${PROJECTION_north} ${OUTTER_north} -X${dx_north} -Y${dy_north} 
        gmt clip -N << END
-80.75 ${lat_in} 
-120 ${lat_in} 
-120 ${lat_out}
-80.75 ${lat_out}
END
            gmt grdimage @earth_relief_01m -Crelief.cpt -I+d
            gmt coast -Dh -Wthinnest,black -Ir/faint,lightblue
        gmt clip -C
    gmt coast -Q

    gmt colorbar -Cice.cpt -DJMR+w${HEIGHT_map}c/0.5c -G0/100 -Np -Bxa10f1g10 -By+l"%"               
    gmt colorbar -Crelief.cpt -DJBC+w${WIDTH_map}c/0.5c+o$(gmt math -Q ${shift} 2 DIV =)c/-1c+h+e -G-4000/1000 -Ba2000f500g1000+l"Topography" -By+l"m" -I                    


    # # # # # SOUTH HEMISPHERE
    ./gmt_plotf.sh ${PROJECTION_south} ${INNER_south} $(gmt math -Q ${shift} ${dx_north} SUB =) $(gmt math -Q ${dy_north} -1 MUL =)

    # Add the piece of Antarctica peninsula that is outside the inner circle
    gmt coast -EAQ+c ${PROJECTION_south} ${OUTTER_south} -X${dx_south} -Y${dy_south} 
        gmt grdimage @earth_relief_01m -Crelief.cpt -I+d
        gmt coast -Dh -Wthinnest,black -Ir/faint,lightblue
    gmt coast -Q

gmt end
EOF

# # # # # # # # # # # # # # # # # # # # # # # # # #
#                  /POST.SH                       #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
#                 SCRIPT.SH\                      #
# # # # # # # # # # # # # # # # # # # # # # # # # #

# Download the data (has to be split due to download size restrictions)
mkdir ./DOWNLOAD
cd DOWNLOAD
for yr in {2013..2016}; 
do
    python download.py $yr
done

# [Extract] the archives and [sort & rename] the files
mkdir ../north
mkdir ../south

# # [extract]
for ff in $(ls -v *zip)
do
    unzip $ff 
done

# # [sort & rename] : north
for nh in $(ls -v *NH*);
do
    mv $nh ../north/$(echo $nh | awk -F[-] '{printf "north_%s.nc",$7}')
done

# # [sort & rename] : south
for sh in $(ls -v *SH*);
do
    mv $sh ../south/$(echo $sh | awk -F[-] '{printf "south_%s.nc",$7}')
done

# GMT TIME !
cd ..  

gmt movie main.sh -Iinclude.sh -Sbpre.sh -Sfpost.sh -Npoles_mov -C1080p -Tlist.txt -D24 -Fmp4 -Pb+w1c+jMC+o0c/3.25c

# Clean the files created "on the fly"
rm -f include.sh pre.sh post.sh main.sh gmt_p*.sh download.py

# # # # # # # # # # # # # # # # # # # # # # # # # #
#                 /SCRIPT.SH                      #
# # # # # # # # # # # # # # # # # # # # # # # # # #