Sunday, March 30, 2008

Script to create sites map from GMT.

#!/bin/sh
#sh_sitemap
# Create site map from sites.defaults files using GMT
#
# Usage: sh_sitemap [ -expt EXPT ] [ -dir ${expt} ] [ -out ofile ]
# [ -check_dc dc.sites.file ]
# default to plot all sites in sites.defaults
#
# Required Files:
# ${expt}/tables/process.defaults
# ${expt}/tables/sites.defaults
# ${expt}/tables/itrfNN*.apr (defined process.defaults)
#
# Input coordiantes are from lfile. or *.apr file.
# lfile. format:
#Epoch 2000.0956: From file itrf00_nafd.apr
#GRAS GRAS GPS N43 33 45.22598 E 6 55 14.05892 6369273.3773 Ref. Epoch 1997.0000 GRAS_GPS
#TOUL TOUL GPS N43 22 7.06826 E 1 28 50.73220 6368233.4576 Ref. Epoch 1997.0000 TOUL_GPS
#7604 7604 GPS N48 13 0.23097 W 4 30 13.78000 6366324.4534 Ref. Epoch 1997.0000 7604_GPS
# apr file format:
#* SYSTEM 1 Coordinates in SYSTEM 2 Frame
#* Name X (m) Y (m) Z (m) Xdot Ydot Zdot Epoch
# GRAS_GPS 4581691.0120 556114.6800 4389360.6960 0.00261 0.02073 -0.00653 1997.0000
# TOUL_GPS 4627846.1280 119629.1780 4372999.7230 0.00232 0.02054 -0.00789 1997.0000
# 7604_GPS 4228877.0780 -333104.1790 4747181.0000 0.00384 0.02027 -0.00663 1997.0000

#Algorithm
# +Get sites names from ${expt}/tables/sites.defaults
# +Get sites coordinates from ${expt}/tables/${itrf}.apr
# ${irtf}.apr filename obtained from ${expt}/tables/process.defaults
# +Plot the sites coordinates with sites names as lables.


#Modification
# + SEP-14-2007 Check the site data availablility from Data Center site list file.
# say, sites.kasi

# initializing parameters
# default to plot all sites coordinates
expt=
# default expt root = current directory
expt_root=`pwd`
# output plot in current directory (hard-wired)
expt_out=`pwd`
# default output name
ofile=sitmap.ps
# default data center site file (to be KASI)
check_dc=site.kasi

while [ "$1" != "" ]; do
case $1 in
-dir)
expt_root=$2
;;
-expt)
expt=$2
;;
-out)
ofile=$2
;;
-check_dc)
check_dc=$2
;;
*)
echo "Invalide option: $1"
exit
;;
esac
shift 2
done

if [ ! -d ${expt_root}/tables ]; then
echo "No ${expt_root}/tables directory exist! Exiting..."
exit
fi


# get sites
file=${expt_root}/tables/sites.defaults
if [ ! -f ${file} ]; then
echo "No ${file} exist! Exiting..."
exit
fi
sites=`grep "^ " ${file} | grep ${expt} | awk '{print substr($1,1,4)}'`
echo "Using Sites: ${sites}"


# site coordiantes files
# use gapr_t_l to convert ${expt}/tables/itrfNN*.apr to old-style lfile.
#


# get apr file
file=${expt_root}/tables/process.defaults
coord_file=`grep aprf ${file}`
# There maybe serveal lines of apr files. Use the last one.
coord_file=`grep aprf ${file}`
coord_file=`echo ${coord_file} | awk '{ for (i=1;i<=NF;i++) {if (i==NF) print $i} }'`
coord_file=`sed 's/#.*//' ${file} | grep aprf | tail -1 | awk '{print $4}'`

#> awk '/aprf/' process.defaults
# set aprf = itrf00_nafd.apr
# set aprf = itrf05.apr
#> sed -n '/line/p' process.defaults
#> sed -n '/aprf/p' process.defaults
# set aprf = itrf00_nafd.apr
# set aprf = itrf05.apr
#> sed -n '/aprf/p' process.defaults | sed -n '$p'
# set aprf = itrf05.apr

coord_file=${expt_root}/tables/${coord_file}
echo ${coord_file}
#exit
# convert .apr to lfile.
#cd ${expt_root}/tables/
#pwd
gapr_to_l ${coord_file} .lf "" 2000.0
coord_file=.lf
#
i=0
tmp_llr=.llr
tmp_lbl=.lbl #labels (site names)
for site in ${sites}; do
echo "Searching coordinates for $site .."
coords=`grep -i ${site}_gps ${coord_file} | awk '{print $1,substr($0,19,14),substr($0,35,15),substr($0,18,1),substr($0,34,1)}' | awk '{ if($9 == "W") {printf("%f", -($5+$6/60+$7/3600))} else {printf("%f", ($5+$6/60+$7/3600))} ;print " "; if($8=="S"){printf("%f",-($2+$3/60+$4/3600))}else{printf("%f",$2+$3/60+$4/3600)} }'`
#echo $coords
#coords_lbl=`grep -i ${site}_gps ${coord_file} | awk '{print $1,substr($0,19,14),substr($0,35,15)}' | awk '{print ($5+$6/60+$7/3600) ,($2+$3/60+$4/3600), "12 0 1 CM ", $1}'` #ERROR
coords_lbl=`grep -i ${site}_gps ${coord_file} | awk '{print $1,substr($0,19,14),substr($0,35,15),substr($0,18,1),substr($0,34,1)}' | awk '{ if($9 == "W") {printf("%f", -($5+$6/60+$7/3600))} else {printf("%f", ($5+$6/60+$7/3600))} ;print " "; if($8=="S"){printf("%f",-($2+$3/60+$4/3600))}else{printf("%f",$2+$3/60+$4/3600)};print " 12 0 1 CM",$1 }'`

if [ "${coords}" != "" ]; then
if [ $i -eq 0 ]; then
echo ${coords} > ${tmp_llr}
echo ${coords_lbl} > ${tmp_lbl}
i=1
else
echo ${coords} >> ${tmp_llr}
echo ${coords_lbl} >> ${tmp_lbl}
fi
else
echo " -Not found."
fi

done
echo "___"
cat ${tmp_llr}
#exit

# plot sites names
#ranges=max/min
rng=`minmax -C ${tmp_llr}`
xmin=`echo ${rng} | awk '{print $1}'`
xmax=`echo ${rng} | awk '{print $2}'`
ymin=`echo ${rng} | awk '{print $3}'`
ymax=`echo ${rng} | awk '{print $4}'`
ix=`echo $xmin $xmax | awk '{print int(int(($2-$1)/10)/5)*5}'`
iy=`echo $ymin $ymax | awk '{print int(int(($2-$1)/10)/5)*5}'`
echo $ix ":" $iy

R=`minmax -I1 ${tmp_llr}`
echo $R

#pscoast $R -I1 -Na -Ggray -JM6i -Wthin -B.5g.5f.5/.5g.25f.25WSen -V -W > ${ofile}
pscoast $R -I1 -Na -Ggray -JM10i -Wthin -B${ix}/${iy}:."SITE MAP":WSEN -V -W -K > ${ofile}

psxy ${tmp_llr} -R -J -O -K -Sc.1i -Gyellow >> ${ofile}

echo "___"
cat ${tmp_lbl}
pstext -R -J -O -N -Gred ${tmp_lbl} >> ${ofile}


\rm -f ${tmp_llr} ${tmp_lbl} .lf fort.*

No comments: