Grillado de superficies

GMT progs : surfer, grdsample
UNIX progs : echo, awk, case, read
En esta sección daremos un ejemplo de grillado de superficie. surfer hace una interpolación forzando a la superficie a pasar por todos los puntos de la grilla, en su proyección (z). El algoritmo que implementa surfer es una versión mejorada al de mínima curvatura. Para más detalles ver Smith and Wessel [1990]
#!/bin/bash

read -p " ingresa mes : " f
case $f in

        enero)  d=0 ;;
        febrero) d=1 ;;
        marzo) d=2 ;;
        abril) d=3 ;;
        mayo) d=4 ;;
        junio) d=5 ;;
        julio) d=6 ;;
        agosto) d=7 ;;
        septiembre) d=8 ;;
        octubre) d=9  ;;
        noviembre) d=10 ;;
        diciembre) d=11 ;;
        *) echo " ERROR!   " 
	echo " ingresa mes con minuscula, por ejemplo  septiembre "
	exit         ;;
esac

more prc.txt |awk '{if( NF == 5) print $0}' | sed s'/\-9.96921e+36/NaN/g' > datos.txt
more datos.txt |grep "time\[$d\]" |sed 's/\=/ /g' |awk '{if( $4 > -57 && $4 < -13 && $6 \
 > 279 && $6 < 315)
print $6-360,$4,$8}' > grilla.xyz


gmtset PLOT_DEGREE_FORMAT=dddF

# definiendo argumentos
	
	region="-81/-45/-58/-13"
	salida="grafo"
	proyeccion="M13"
	topografia="/vox/201001/linux/topografia/topo_6.2.img"
	marco="a9f9/a8f8NWse"
	illaz="190"
	bordecos="0.01p/0/0/0"
	limfron="-N1/0.5p/0"
	psfile="${salida}.ps"
	paleta="/vox/201001/linux/certamen/GMT_ocean.cpt"


# ajustes iniciales

 img2grd ${topografia} -G${salida}.grd -R${region} -m2.0 -T1 -N1 -V
 psbasemap -B${marco} -J${proyeccion} -R${region} -X3 -Y5 -P -V -K > ${psfile}
 grdgradient ${salida}.grd -G${salida}.int  -A${illaz} -Nt -M
 
# interpolacion

 surface ${grilla}.xyz -G${grilla}.grd -R${region} -I1m -T0.2 -C1 -V
 makecpt -C/vox/201001/linux/certamen/GMT_meanprecip.cpt   -T0/10/1 -I  > ${grilla}.cpt
 grdsample ${salida}.int -Gtopo.int -R${region} -I1m -V -F
 grdsample ${grilla}.grd -Gmes2.grd -R${region} -I1m -V -F
 grdimage mes2.grd -C${grilla}.cpt -B${marco} -J${proyeccion} -Itopo.int -R${region}\
  -V -O -K >> ${psfile}
 grdcontour mes2.grd -B${marco} -J${proyeccion} -C50 -A1+k0/0/0+g255/255/255 -R${region} \
 -W0.5p/0/0/0 -Gd5c -O -K >> ${psfile}

# ajustes finales
 
 pscoast -B${marco} -J${proyeccion} -R${region}  -P -V -Df -S0 -W${bordecos}\
 ${limfron} -O -K >> ${psfile}
 echo 49.7W 49.7S 1 2 | psxy -R${region} -J${proyeccion} -O -K -Sr -G255\
  -W0.5  >> ${psfile}
 echo 50.5W 50S | mapproject -R279/315/-58/-13 -JM13  > co
 pos=` awk '{printf "%sc/%sc\n", $1, $2}' co `
 max="a10f1/a10f1:mm:"
 psscale  -D$pos/4/0.3 -C${grilla}.cpt ${verbose} -B${max} -I -O  >> ${psfile}

rm *.grd *.int *.xyz  co  mes.cpt
Explicación
1.
El archivo a procesar prc.txt es extraido de un netCDF, se suguiere revisar detalladamente en que consiste el .txt. Este .txt contiene la precipitación media para cada mes del año entre 1979 - 2009. La gráfica será sobre la región de America del Sur
2.
En el script la combinación de read y case dan la eleción del mes a graficar
3.
Una vez que tenemos nuestra variable (mes = $d) hacemos una depuración de los datos a fin de poder procesar. Para esto utilizamos los comandos sed, awk, more
4.
Definimos argumentos en GMT, además de hacer una subgrilla con el comando img2grd
5.
surface interpola una grilla a partir de un archivo de coordenadas (x,y,z), donde los parámetros -T, -I, -C dan el factor de tensión de la superficie, el espaciado de grilla y el nivel de convergencia de las iteraciones ligadas a la interpolación, respectivamente
6.
Luego creamos una nueva paleta a partir de una madre asignando valores extremos deacuerdo a nuestro set de datos xyz
7.
grdsample se usa para re muestrear una grilla , creando otra con diferente espaciemiento u otra subregión. Esto se hace para el archivo .grd asi como también para el archivo que contiene la iluminación
8.
Aplicamos contornos a gusto
9.
Finalmente graficamos el borde de costa y usando nustros avances de UNIX hacemos unos trucos para modificar la paleta de colores que irá al lado del grafico

gonzalo 2011-09-16