GMT Scatterplot

This script can be used to draw a mediawiki:Scatterplot using data provided in a comma seperated text file.

Note

You need the following to run a GMT program successfully

• Install GMT
• A UNIX like environment with standard programs like awk and sed. If you use Linux you are already set. If you use windows Install Cygwin -- it's very simple to do this.

You need the following to use the Postscript output of GMT

• Free utility to view Postscript files -- Ghostview or GSview.

Script

filename: x-y.bash

#!/bin/bash
# (c) 2006, Assela Pathirana under GNU GPL 2.0 (http://www.gnu.org/copyleft/gpl.html)
#############
#Do some input checking and help out the first time users
sourcef="$1.txt"
if [[ $# -lt 1 ]] 
then
  echo 'to create a scattergram from the datafile named <arg1>''.txt'' and to save it in <arg1>.eps'
  echo "The data file <arg1>.txt should  have the data in two columns, seperated by commas. The first line of input should be the axes titles, not data. "
  echo 'first run with only \<arg1> as file name. '
  echo "Then after seeing the range of data to plot - run with all five arguments <arg1> <dx> <dy>"
  exit 1
else
  if [[ ! -f "$sourcef" ]]; then
    echo "There is no file named $sourcef (did you forget to give filename witout .txt part?) " 
    exit 1
  fi
fi
if [ $# -ge 3 ]
then
  ytitle=`head -n1 $sourcef |awk 'BEGIN{FS=","}{print $2}'`
  xtitle=`head -n1 $sourcef |awk 'BEGIN{FS=","}{print $1}'`
	inilt=-I$2/$3
	b="a$2f$2g$2:$xtitle:/a$3f$3g$3:$ytitle:WSne"
  echo "using -B option: "$b
else
  if [ $# -eq 1 ] 
  then
    inilt=-I1E-15/1E-15
    b="a1f1g1WSne"
  else
#now we are in a strange position. 
    echo "Either give one argument or all five! I quit. "
    exit 1
  fi
fi
#Simple check for gmt commands
if [[ ! -f  "$GMTHOME/bin/gmtset" ]] ; then 
 echo "no GMT environment present. Install GMT and set paths properly before running this script"
 exit 1
fi

#We are In business
#clear any garbage from prvious gmt runs/attempts
rm -f .gmtcommands* .gmtdefaults*
#set some constants
gmtset LABEL_FONT_SIZE 18
outf="$1.eps"
trcol="1pt/255/0/0"
bicol="1pt3_3:0/0/0/255"
symcol="255/0/0"
symsym="C.4c"
trendfile="$1.tr.txt"
# calculate a linear(N2) fit for the datafile who has one header rows. Write only 
trend1d $sourcef   -N2 -Fxmy -H1 |sort -bg > $trendfile

# calculate the coefficient of determination
ybar=`cat $trendfile |awk '{ya=ya+$3;ct++}END{print ya/ct}'`
CD=`cat $trendfile |awk -v ybar=$ybar '{ymb=ymb+($3-ybar)^2;ymh=ymh+($3-$2)^2;ct++}END{printf "%4.2f", (ymb-ymh)/ymb}'`
#find a plotting range
minmax $sourcef -H1 $inilt -C > $$
xmin=`awk  -v dd=$2 '{print $1-dd}' $$`
ymin=`awk  -v dd=$3 '{print $3-dd}' $$`
xmax=`awk  -v dd=$2 '{print $2+dd}' $$`
ymax=`awk  -v dd=$3 '{print $4+dd}' $$`
\rm -f $$
range='-R'$xmin'/'$xmax'/'$ymin'/'$ymax'/'
echo "Using the range : $range"
echo "now plotting ..."
psxy -H1 $sourcef $range -W.5p/0/0/0 -G$symcol -S$symsym  -JX10c -B"$b" -P  -K > $outf
psxy -H1 $trendfile $range -W1.5p/$trcol   -JX10c -B"$b" -P -O -K >> $outf

#legend
psxy -R0/1/0/1 -JX -O -K -G240 -L -W.75p << END >> $outf 
.055 .955
.405 	.955 
.405	.705 
.055	.705 
END
echo "0.09 0.90" | psxy -R -JX -O -K -S$symsym -G$symcol -W.5p/0/0/0  >> $outf 
echo "0.16 0.90 14 0.0 1 5 Values" | pstext -R -JX -O -K >> $outf 
psxy -R -JX -O -K -W$trcol -G0 <<END >>$outf
0.07	0.84
0.14 	0.84
END
echo "0.16 0.84 14 0.0 1 5 Trend" | pstext -R -JX -O -K >> $outf 

psxy -R0/1/0/1 -JX -O -K -G240 -L -W.75p << END >> $outf 
.96 .06	
.96 .21 
.61 .21	
.61 .06	
END
echo "0.65 0.135 14 0.0 1 5 C.D.= $CD" | pstext -R -JX -O   >> $outf
#echo "0.16 0.78 14 0.0 1 5 `head -n1 $sourcef|awk 'BEGIN{FS=","};{print $1}'` " | pstext -R -JX -O >> $outf 
#clean up
\rm -f .gmt* $$ $$.* 

Test Data

filename: sample.data.txt

x data title , y data title
-10.681   ,  -4.14932696
-8.292  ,  -0.665670484
-25.482   ,  -3.775538869
-24.293   ,  -0.315713513
-36.929   ,  -4.215689714
-42.547   ,  -2.864580275
-49.196   ,  -8.841218073
-54.941   ,  -4.88540979
-63.0095  ,  -10.46740317
-69.0207  ,  -8.121834461

How to Run

In a UNIX environment (which includes Cygwin and Linux), the user need to have executable permissions on a file to run it. (Here is a UNIX filesystem permissoins tutorial.)

---