GMT Scatterplot

From assela Pathirana
Jump to navigationJump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
Magick .jpg

This script can be used to draw a 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

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.)