#!/bin/bash
################################################################################
# Script que pinta con GMT mapas para José Miguel Viñas.
#
# Marcos Molina Cano <marcosmolina.tiempo@gmail.com>                  26/01/2018
################################################################################
# Para buscar errores.
set -x

# Definición de parámetros de la pasada.
dia="16"
mes="02"
ano="2018"
fecha=${ano}${mes}${dia}
pasada="00"
mesLetra=$(date -u +%b --date="${ano}-${mes}-${dia} ${pasada}:00:00 UTC")

proyInicial=24
proyFinal=168
dt=12

# Directorio donde están los gribs, los NetCDF y las imágenes.
GRIB_DIR=/data/ECMWF/HRES-WAM/${fecha}${pasada}
NC_DIR=${GRIB_DIR}_NC
IMAGES_DIR=${GRIB_DIR}_IMAGES
# Nos aseguramos que existen los directorios.
mkdir -p ${GRIB_DIR} ${NC_DIR} ${IMAGES_DIR}


# Comienzo del trabajo para cada proyección.
for PROY in $(seq --format %03g ${proyInicial} ${dt} ${proyFinal}) ; do

  # Fechas de validez de cada proyección.
  horaVal=$(date -u +%H --date="${ano}-${mes}-${dia} ${pasada}:00:00 UTC ${PROY} hours")
  diaVal=$(date -u +%d --date="${ano}-${mes}-${dia} ${pasada}:00:00 UTC ${PROY} hours")
  mesVal=$(date -u +%m --date="${ano}-${mes}-${dia} ${pasada}:00:00 UTC ${PROY} hours")
  mesValLetra=$(date -u +%b --date="${ano}-${mes}-${dia} ${pasada}:00:00 UTC ${PROY} hours")
  anoVal=$(date -u +%Y --date="${ano}-${mes}-${dia} ${pasada}:00:00 UTC ${PROY} hours")


  # Para el modelo marítimo HRES-WAM ya están en una malla
  # de tipo regular_ll y como todos los campos pertenecen
  # al mismo nivel podemos pasar directamente a NetCDF.
  root=A1P${mes}${dia}${pasada}00${mesVal}${diaVal}${horaVal}001
  gribFile=${GRIB_DIR}/${root}
  ncFile=${NC_DIR}/${root}.nc
  testFile=${gribFile}.ok
  if test ! -e ${testFile} ; then
    if test -s ${gribFile} ; then
      grib_to_netcdf ${gribFile} -o ${ncFile}
      touch ${testFile}
    else
      echo "No se encuentra ${gribFile}. ¿Quieres copiarlo? [S/n]"
#      read respuesta
      respuesta="S"
      if test "${respuesta}" = "" -o "${respuesta}" = "S"   ; then
        echo "Vamos a copiarlo."
        sshpass -p '' scp ecmwf@dali.ogimet.com:/home/ecmwf/${root} ${GRIB_DIR}
#        sshpass -p '' scp ecmwf@dali.ogimet.com:/home/ecmwf/ecmwf_${fecha}${pasada}/${root} ${GRIB_DIR}
      elif test "${respuesta}" = "n" ; then
        echo "Nos vemos."
        exit
      else
        echo "Por favor, responde con la letra S o n."
        exit
      fi
    fi

    # Comprobamos si se ha copiado correctamente el grib.
    if test -s ${gribFile} ; then
      grib_to_netcdf ${gribFile} -o ${ncFile}
      touch ${testFile}
    else
      echo "Problemas al copiar ${gribFile}."
      exit
    fi
  fi

  # Ficheros NetCDF con los datos a pintar y sus correspondientes arreglos.
  #------------------------------------------------------------------------

  # hmax = Altura máxima del oleaje. (m)
  ncFileHmax=${NC_DIR}/hmax.nc
  gmt grdclip ${ncFile}?hmax -G${ncFileHmax} -SrNaN/0

  # mwdU = Componente U de la dirección media del oleaje. (º)
  ncFileMwdu=${NC_DIR}/mwdu.nc
  gmt grdmath 270 ${ncFile}?mwd SUB COSD = ${ncFileMwdu}

  # mwdV = Componente V de la dirección media del oleaje. (º)
  ncFileMwdv=${NC_DIR}/mwdv.nc
  gmt grdmath 270 ${ncFile}?mwd SUB SIND = ${ncFileMwdv}

  # mdw = Dirección media del oleaje en polares planas.
  ncFileMwd=${NC_DIR}/mwd.nc
  gmt grdmath 270 ${ncFile}?mwd SUB = ${ncFileMwd}

  # ci = Sea ice cover (%)
  ncFileCi=/data/ECMWF/ci.nc
  if test ! -s ${ncFileCi} ; then
    echo "No se encuentra ${ncFileCi}"
    exit
  fi

  ###########################################
  #  Comienzo del pintado del mapa con GMT. #
  ###########################################
  # Nombre de la imagen.
  psFile=${GRIB_DIR}_IMAGES/${fecha}${pasada}_${PROY}_hmax_mwd.ps
  pngFile=${GRIB_DIR}_IMAGES/${fecha}${pasada}_${PROY}_hmax_mwd.png

  # Definición de parámetros de la proyección para la opción -J.
  lon0=-65
  lat0=-59
  altitude=1200 # Altura en km.
  azimuth=0  # azimuth is measured to the east of north of view
  tilt=0 # tilt is the upward tilt of the plane of projection
  twist=0 # Further, specify the clockwise twist, Width, and Height of the viewpoint in degrees.
  Width=90 # Anchura en grados.
  Height=${Width} # Altura en grados.
  width=16c # ancho del mapa en centímetros.
  J=G${lon0}/${lat0}/${altitude}/${azimuth}/${tilt}/${twist}/${Width}/${Height}/${width}

  # Definición del área geográfica a pintar.
  R=d

  # Paletas de colores.
  Chmax=/home/marcos/wrf_model/share/cpt/hmax.cpt
  Cci=/home/marcos/wrf_model/share/cpt/ci.cpt

  # Localización y configuración de escalas.
  Dhmax=n0.05/-0.05+w6/0.25+h+ef0.5
  Dci=n0.55/-0.05+w6/0.25+h

  # Estilo de las líneas que delimitan paralelos y meridianos.
  gmt set MAP_GRID_PEN_PRIMARY="0.5p,black,dotted"

  # Altura máxima del oleaje.
  gmt grdimage   -R${R} -J${J} ${ncFileHmax}  -C${Chmax}   -E300 -nb -Q -K    >  ${psFile}

  # Cobertura de hielo marítimo.
  gmt grdimage   -R     -J     ${ncFileCi}?ci -C${Cci}     -E300 -nb -Q -K -O >> ${psFile}

  # Líneas de costa y colorado de continentes.
  gmt pscoast   -R -J -Dh -N1 -N3 -Gdarkgray -A10 -W -K -O --MAP_FRAME_TYPE="plain" >> ${psFile}

  # Escalas.
  gmt psscale   -R -J -D${Dhmax}   -C${Chmax} -S -L -O -K >> ${psFile}
  gmt psscale   -R -J -D${Dci} -C${Cci} -S -L -O -K >> ${psFile}

  # Rótulos superiores.
  echo "Altura m@\341xima y direcci@\363n media del oleaje." \
      | gmt pstext -D-8c/8.8c -N -R -J -F+f12p,Helvetica-Bold,black+cMC+jLB    -K -O -V >> ${psFile}

  echo "ECMWF HRES-WAM (0.125@\260)" \
      | gmt pstext -D-8c/8.2c -N -R -J -F+f12p,Helvetica-Bold,darkgray+cMC+jLB -K -O -V >> ${psFile}

  echo "Inicio: ${dia} ${mesLetra} ${ano}, ${pasada} UTC" \
      | gmt pstext -D+8c/8.8c -N -R -J -F+f12p,Helvetica,darkred+cMC+jRB       -K -O -V >> ${psFile}

  echo "V@\341lido: ${diaVal} ${mesValLetra} ${anoVal}, ${horaVal} UTC (H+${PROY})" \
      | gmt pstext -D+8c/8.2c -N -R -J -F+f12p,Helvetica-Bold,darkred+cMC+jRB  -K -O -V >> ${psFile}

  # Rótulos de las escalas.
  echo "Altura m@\341xima del oleaje (m)"    | gmt pstext -D-4c/-8.2c  -N -R -J -F+f12p,Helvetica,black+cMC -O -K -V >> ${psFile}
  echo "Cobertura de hielo mar@\355timo (%)" | gmt pstext -D+4c/-8.2c  -N -R -J -F+f12p,Helvetica,black+cMC -O -K -V >> ${psFile}

  # Dirección del oleaje marítimo.
#  gmt grdvector -R -J ${ncFileMwdu} ${ncFileMwdv} -Ix8/4 -Sl50 -Q0.2c+a15+e+jc -Wdarkbrown -Gdarkbrown -K -O -Vl --MAP_VECTOR_SHAPE="0"  >> ${psFile}
  gmt grdvector -R -J ${ncFileHmax} ${ncFileMwd} -A -Ix6/3 -Si5 -Q0.15c+a25+e++jc+n30 -Wdarkbrown -Gdarkbrown --MAP_VECTOR_SHAPE="0" -K -O -Vl >> ${psFile}
  #gmt grdvector -R -J ${ncFile}?hmax ${ncFile}?mwd -A -Ix9 -Sl100 -Q0.2c+a15+b+jc -Wblue -Gblue -K -O -Vl --MAP_VECTOR_SHAPE="0" -t50 >> ${psFile}

  # Logos de tiempo.com y ECMWF.
  gmt psimage ~/wrf_model/scripts/GMT/ES.eps    -J -R -Dx11.9c/0.1c+w4c -F+gwhite -K -O -Vl -t20 >> ${psFile}
  gmt psimage ~/wrf_model/scripts/GMT/ecmwf.eps -J -R -Dx0.1c/0.1c+w4c  -F+gwhite -K -O -Vl -t20 >> ${psFile}

  # Se añade marco.
  gmt psbasemap -R -J -Bp5 -Bg5 -O -K --MAP_FRAME_TYPE="inside" --MAP_FRAME_AXES="WN" --MAP_ANNOT_OFFSET="2p" --FONT_ANNOT_PRIMARY="8p,Helvetica-Bold,black" >> ${psFile}

  # Se machagan logos de tiempo.com y ECMWF sin transparencia.
  gmt psimage ~/wrf_model/scripts/GMT/ecmwf.eps -J -R -Dx0.1c/0.1c+w4c            -K -O -Vl      >> ${psFile}
  gmt psimage ~/wrf_model/scripts/GMT/ES.eps    -J -R -Dx11.9c/0.1c+w4c              -O -Vl      >> ${psFile}

  # Se convierte a png.
  gmt psconvert -Z -A0.2c/0.1c -P -Tg -Qg4 -Qt4 ${psFile}

  # Se transforma.
  convert ${pngFile} -resize 800 ${pngFile}

  # pdf2svg ${WORK_DIR}/genper.pdf ${WORK_DIR}/genper.svg

  # Limpiamos.
  rm -f ${ncFileHmax} ${ncFileMwdu} ${ncFileMwdv}

done