[cig-commits] [commit] master: add the GMT script grdview.sh to ease the creation of 3D perspectives. (294099a)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Jun 13 22:32:35 PDT 2014
Repository : https://github.com/geodynamics/relax
On branch : master
Link : https://github.com/geodynamics/relax/compare/a5df9b4d88e7dcc311b642c553d528f43ecf2df0...672813afe9de9ae01c5c16181fed308f2b5203b2
>---------------------------------------------------------------
commit 294099a945f957114e748034c508fb326c713e9d
Author: Sylvain Barbot <sbarbot at ntu.edu.sg>
Date: Fri Jun 13 22:31:21 2014 -0700
add the GMT script grdview.sh to ease the creation of 3D perspectives.
>---------------------------------------------------------------
294099a945f957114e748034c508fb326c713e9d
util/{grdmap.sh => grdview.sh} | 186 ++++++++++++++++++-----------------------
1 file changed, 81 insertions(+), 105 deletions(-)
diff --git a/util/grdmap.sh b/util/grdview.sh
similarity index 61%
copy from util/grdmap.sh
copy to util/grdview.sh
index 95a3d68..187bd24 100755
--- a/util/grdmap.sh
+++ b/util/grdview.sh
@@ -1,11 +1,11 @@
#!/bin/bash
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
-# script grdmap.sh
-# display a map view of GMT .grd files with
+# script grdview.sh
+# display a view of GMT .grd files with
# overlay from user-provided GMT scripts.
#
-# Sylvain Barbot (01/01/07) - original form
+# Sylvain Barbot (06/13/07) - original form
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
set -e
@@ -15,28 +15,30 @@ cmdline=$*
trap 'echo $self: Some errors occurred. Exiting.; exit' ERR
usage(){
- echo "usage: $self [-b -5/5/-5/5] [-p -1.5/1.5/0.1] [...] file1.grd ... fileN.grd"
+ echo "usage: $self [-b -5/5/-5/5/-5/5] [-E 180/90] [-p -1.5/1.5/0.1] [...] file1.grd ... fileN.grd"
echo "or"
- echo " $self [-b -5/5/-5/5] [-p -1.5/1.5/0.1] [...] dir1/index1 ... dirN/indexN"
+ echo " $self [-b -5/5/-5/5/-5/5] [-E 180/90] [-p -1.5/1.5/0.1] [...] dir1/index1 ... dirN/indexN"
echo ""
echo "options:"
echo " -b bds sets the map bound to bds"
echo " -c palette_name [default my_jet]"
- echo " -e file.sh runs file.sh to add content to map"
+ echo " -e efile.sh runs efile.sh to add content to map"
echo " -g switch to geographic projection (longitude/latitude)"
- echo " -h display this error message and exit"
+ echo " -h display this help message and exit"
echo " -i file.grd illuminates the grd image with file.grd"
echo " -o output directory"
echo " -p palette sets the color scale bounds to palette"
echo " -r reverse the y-axis"
echo " -s step sets the distance between vectors"
echo " -u defines the color scale unit"
- echo " -v vector sets the vector scale to vector"
+ echo " -v vertical exxageration"
echo " -t tick interval"
echo " -T title header"
echo " -x do not display map (only create .ps file)"
- echo " -C interval plots contours every interval distance"
- echo " -E file.ps only plot base map with extra scripts"
+ echo " -E Sets the viewpoint's azimuth and elevation for perspective view [180/90]"
+ echo " -G drapefile"
+ echo " -N base level"
+ echo " -O file.ps only plot base map with extra scripts"
echo " -J overwrites the geographic projections (for -J o the -b option is relative)"
echo " -Y shift the plot vertically on the page"
echo ""
@@ -50,21 +52,19 @@ usage(){
my_gmt(){
if [ -e "$U3" ]; then
- grdimage $U3 -R$bds -J${PROJ} -Q \
- $AXIS \
- -K -C$colors -P -X1.2i -Y${YSHIFT}i $illumination \
+ grdview $U3 -R$bds $PROJ $PROJZ -Qc $Eset $Nset \
+ $AXIS $Gset \
+ -K -C$colors -P -X1.2i -Y${YSHIFT}i $iset \
> $PSFILE
- if [ "$Cset" == "-C" ]; then
- grdcontour -K -O $U3 -W0.3p/100/100/100 -R$bds -J${PROJ} -S -C$contour -P >> $PSFILE
- fi
else
- psbasemap -R$bds -J${PROJ} \
- $AXIS \
+ echo $self: psbasemap -R$bds $PROJ $PROJZ $AXIS $Eset -K -P
+ psbasemap -R$bds $PROJ $PROJZ \
+ ${AXIS}wsne $Eset \
-K -P -X1.2i -Y${YSHIFT}i \
> $PSFILE
fi
- # running all required subprograms
+ # add extra layers
for subprog in $EXTRA; do
if [ "" != "$U3" ]; then
OPTIONP="-p $U3"
@@ -73,49 +73,18 @@ my_gmt(){
fi
if [ -e "$selfdir/$subprog" ]; then
#echo $self: running $subprog $PSFILE $bds $VECTOR $U3 $HEIGHT
- eval "$subprog $gset -b $bds -v $VECTOR $OPTIONP -H $HEIGHT $PSFILE"
+ eval $subprog $gset -b $bds $OPTIONP -H $HEIGHT -E ""$Eset"" -J ""$PROJ"" $PSFILE
else
if [ -e "$subprog" ]; then
- #echo $self: running $subprog $PSFILE $bds $VECTOR $U3 $HEIGHT
- eval "$subprog $gset -b $bds -v $VECTOR $OPTIONP -H $HEIGHT $PSFILE"
+ EsetP='\"$Eset\"'
+ PROJP='\"$PROJ\"'
+ PROJZP='\"$PROJZ\"'
+ #echo $self: running $subprog $gset -b $bds $OPTIONP -H $HEIGHT -E $EsetP -J $PROJP -Z $PROJZP $PSFILE
+ eval $subprog $gset -b $bds $OPTIONP -H $HEIGHT -E $EsetP -J $PROJP -Z $PROJZP $PSFILE
fi
fi
done
- # plotting vectors
- if [ -e "$U1" ]; then
-
- #echo $self": found "$U1": plotting vectors"
- echo $self": Using VECTOR="$VECTOR", STEP="$ADX
-
- # arrowwidth/headlength/headwidth
- # grdvector crashes with wrong sampling
- # -Q0.51c/0.8c/0.7cn1.0c \
- # -Q0.3c/0.5c/0.4cn1.0c \
- grdvector $U1 $U2 -K -J${PROJ} -O -R$bds -I$ADX/$ADX \
- -Q0.05/0.2/0.08n0.3c \
- -S$VECTOR -W0.5p/0/0/0 \
- -G255/255/255 \
- >> $PSFILE
- fi
-
- # plot the vector legend if $VECTOR is set
- if [ "$VECTOR" != "-" ]; then
- UL=`echo $bds | awk -F "/" '{print $1,$4}' `
- pstext -O -K -J${PROJ} -N -R$bds \
- -G0/0/0 -Yr0.3i \
- <<EOF >> $PSFILE
-$UL 14 0 4 LM $SIZE $unit
-EOF
- psxy -O -K -J${PROJ} -R$bds -N \
- -W0.5p/0/0/0 -Xr0.9i \
- -Sv0.05/0.2/0.08n0.3c \
- <<EOF >> $PSFILE
-$UL 0 1
-EOF
- REVERT="-Xr-0.9i -Yr-0.3i"
- fi
-
if [ -e "$colors" ]; then
#-Q0.20c/1.0c/0.4cn1.0c \
psscale -O -K -B$PSSCALE/:$unit: -D3.5i/-0.8i/7.1i/0.2ih \
@@ -125,7 +94,10 @@ EOF
rm -f $colors
fi
- echo 0 0 | psxy -O -J${PROJ} -R$bds -Sp0.001c >> $PSFILE
+ psbasemap -R$bds $PROJ $PROJZ \
+ ${AXIS}WSNEZ+ $Eset -O -P >> $PSFILE
+
+ #echo 0 0 | psxy -O $PROJ -R$bds -Sp0.001c >> $PSFILE
}
@@ -140,31 +112,32 @@ gmtset PAPER_MEDIA archA
libdir=$(dirname $0)/../share
EXTRA=""
-while getopts "b:c:e:ghi:o:p:v:s:t:T:u:xrC:E:J:Y:" flag
-do
+while getopts "b:c:e:ghi:o:p:v:s:t:T:u:xrE:G:N:O:J:Y:" flag; do
case "$flag" in
b) bset=1;bds=$OPTARG;;
c) cset="-c";carg=$OPTARG;;
e) eset=1;EXTRA="$EXTRA $OPTARG";;
g) gset="-g";;
h) hset=1;;
- i) iset=1;illumination="-I$OPTARG";;
+ i) iset="-I$OPTARG";;
o) oset=1;ODIR=$OPTARG;;
p) pset=1;P=$OPTARG;PSSCALE=`echo $OPTARG | awk -F"/" 'function abs(x){return x<0?-x:x}{print abs($2-$1)/6}'`;;
r) rset=1;;
s) sset=1;ADX=$OPTARG;;
- t) tset=1;tick=$OPTARG;;
+ t) tset=1;tick=$OPTARG;tickz=$OPTARG;;
T) Tset=1;title=$OPTARG;;
u) uset=1;unit=$OPTARG;;
- v) vset=1;SIZE=$OPTARG;VECTOR=$OPTARG"c";;
+ v) vset=1;verticalexxageration=$OPTARG;;
x) xset=1;;
- C) Cset="-C";contour=$OPTARG;;
- E) Eset=1;PSFILE=$(dirname $OPTARG)/$(basename $OPTARG .ps).ps;;
- J) Jset=1;PROJ=$OPTARG;;
+ E) Eset="-E$OPTARG";;
+ G) Gset="-G$OPTARG";;
+ N) Nset="-N$OPTARG";;
+ O) Oset=1;PSFILE=$(dirname $OPTARG)/$(basename $OPTARG .ps).ps;;
+ J) Jset=1;PROJ="-J$OPTARG";;
Y) Yset=1;Yshift=$OPTARG;;
esac
done
-for item in $bset $cset $iset $oset $pset $vset $sset $tset $Tset $uset $Yset $Cset $Eset $Jset $EXTRA;do
+for item in $bset $cset $iset $oset $pset $sset $tset $vset $Tset $uset $Yset $Eset $Gset $Nset $Oset $Jset $EXTRA;do
shift;shift
done
for item in $gset $hset $xset $rset;do
@@ -177,7 +150,7 @@ done
if [ "$vset" != "1" ]; then
# unused value but preserve the number of elements in call to subroutine
- VECTOR="-"
+ verticalexxageration=1
fi
# vertical shift of plot
@@ -198,8 +171,13 @@ else
fi
fi
+# azimuth and elevation
+if [ "$Eset" == "" ]; then
+ Eset="-E180/90"
+fi
+
# usage
-if [ $# -lt "1" -a "$Eset" != "1" ] || [ "$hset" == "1" ] ; then
+if [ $# -lt "1" -a "$Oset" != "1" ] || [ "$hset" == "1" ] ; then
usage
else
echo $self: using colorfile $cptfile
@@ -207,15 +185,13 @@ fi
# loop over grd files
-while [ "$#" != "0" -o "$Eset" == "1" ];do
+while [ "$#" != "0" -o "$Oset" == "1" ];do
if [ "$1" != "" ]; then
WDIR=`dirname $1`
INDEX=`basename $1`
# default names
- U1=$WDIR/$INDEX-east.grd
- U2=$WDIR/$INDEX-north.grd
U3=$WDIR/$INDEX-up.grd
if [ -e "$U3" ]; then
@@ -224,14 +200,14 @@ while [ "$#" != "0" -o "$Eset" == "1" ];do
U3=$WDIR/$INDEX
fi
if [ "$bset" != "1" ]; then
- bds=`grdinfo -C $U3 | awk '{s=1;print $2/s"/"$3/s"/"$4/s"/"$5/s}'`
+ bds=`grdinfo -C $U3 | awk '{s=1;print $2/s"/"$3/s"/"$4/s"/"$5/s"/"$6/s"/"$7/s}'`
fi
# tick marks
if [ "$tset" != "1" ]; then
tick=`echo $bds | awk -F "/" '{s=1;print ($2-$1)/s/4}'`
fi
- tickf=`echo $tick | awk '{print $1/2}'`
+ tickf=`echo $tick | awk '{print $1/4}'`
echo $self": Using directory "$WDIR", plotting index "$INDEX
@@ -265,37 +241,32 @@ while [ "$#" != "0" -o "$Eset" == "1" ];do
else
grd2cpt $U3 -C$cptfile -Z -T= > $colors
fi
- fi
-
- if [ -e "$U1" ]; then
- if [ "$vset" != "1" ]; then
- MAX1=`grdinfo $U1 -C | awk '{if (sqrt($7^2)>sqrt($6^2)){print sqrt($7^2)}else{print sqrt($6^2)}}'`
- MAX2=`grdinfo $U2 -C | awk '{if (sqrt($7^2)>sqrt($6^2)){print sqrt($7^2)}else{print sqrt($6^2)}}'`
- SIZE=`echo "$MAX1 $MAX2"| awk '{if (0==$1 && 0==$2){print 1}else{print ($1+$2)*0.95}}'`
- VECTOR=$SIZE"c"
- fi
- if [ "$sset" != "1" ]; then
- ADX=`grdinfo $U2 -C | awk '{printf "5*%11.9f\n", $8}' | bc`
- fi
+ else
+ TRIANGLES=`grdinfo -C $U3 | \
+ awk -v m=$m 'function max(x,y){return (x>y)?x:y}; \
+ function abs(x){return (0<x)?x:-x}{if (m<max(abs($6),abs($7))){print "-E"}}'`
fi
if [ "$gset" != "-g" ]; then
if [ "$Jset" == "" ]; then
# Cartesian coordinates
- HEIGHT=`echo $bds | awk -F "/" '{printf("%fi\n",($4-$3)/($2-$1)*7)}'`
+ HEIGHT=`echo $bds | awk -F "/" '{printf("%fi\n",($4-$3)/($2-$1)*5)}'`
+ SCALE=`echo $bds | awk -F "/" -v s=$verticalexxageration '{printf("%f\n",5/($2-$1)*s)}'`
if [ "$rset" != "1" ]; then
- PROJ="X7i/"$HEIGHT
+ PROJ="-JX5i/"$HEIGHT
+ PROJZ=-Jz$SCALE
else
- PROJ="X7i/"-$HEIGHT
+ PROJ="-JX5i/"-$HEIGHT
+ PROJZ=-Jz$SCALE
fi
else
HEIGHT="-"
fi
- AXIS=-Bf${tickf}a${tick}:"":/f${tickf}a${tick}:""::."$title":WSne
+ AXIS=-Bf${tickf}a${tick}:"":/f${tickf}a${tick}:"":/f25a100:""::."$title":
else
# geographic coordinates
- HEIGHT=7i
- PROJ="M$HEIGHT"
+ HEIGHT=5i
+ PROJ="-JM$HEIGHT"
AXIS=-B${tick}:"":/${tick}:""::."$title":WSne
fi
@@ -315,30 +286,35 @@ while [ "$#" != "0" -o "$Eset" == "1" ];do
colors=$WDIR/palette.cpt
fi
if [ "$bset" != "1" ]; then
- echo "$self: -b option should be set with the -E option. exiting."
+ echo "$self: -b option should be set with the -O option. exiting."
echo ""
usage
else
# tick marks
if [ "$tset" != "1" ]; then
- tick=`echo $bds | awk -F "/" '{s=1;print ($2-$1)/s/7}'`
+ tick=`echo $bds | awk -F "/" '{s=1;print ($2-$1)/s/5}'`
+ tickz=`echo $bds | awk -F "/" '{s=1;print ($6-$5)/s/5}'`
fi
- tickf=`echo $tick | awk '{print $1/2}'`
+ tickf=`echo $tick | awk '{print $1/4}'`
+ tickzf=`echo $tickz | awk '{print $1/4}'`
# Cartesian vs geographic coordinates
if [ "$gset" != "-g" ]; then
- HEIGHT=`echo $bds | awk -F "/" '{printf("%fi\n",($4-$3)/($2-$1)*7)}'`
+ HEIGHT=`echo $bds | awk -F "/" '{printf("%fi\n",($4-$3)/($2-$1)*5)}'`
+ SCALE=`echo $bds | awk -F "/" -v s=$verticalexxageration '{printf("%f\n",5/($2-$1)*s)}'`
if [ "$rset" != "1" ]; then
- PROJ="X7i/"$HEIGHT
+ PROJ=-JX5i/$HEIGHT
+ PROJZ=-Jz$SCALE
else
- PROJ="X7i/"-$HEIGHT
+ PROJ=-JX5i/-$HEIGHT
+ PROJZ=-Jz$SCALE
fi
- AXIS=-Ba${tick}:"":/a${tick}:""::."$TITLE":WSne
+ AXIS=-Ba${tick}:"":/a${tick}:"":/f${tickzf}a${tickz}:""::."$TITLE":
else
- HEIGHT=7i
- PROJ="M0/0/$HEIGHT"
- AXIS=-B${tick}:"":/${tick}:""::."$TITLE":WSne
+ HEIGHT=5i
+ PROJ="-JM0/0/$HEIGHT"
+ AXIS=-B${tick}:"":/${tick}:"":/f50:""::."$TITLE":WSne
fi
fi
fi
@@ -346,8 +322,8 @@ while [ "$#" != "0" -o "$Eset" == "1" ];do
# color bounds
if [ "$pset" == "1" ]; then
makecpt -C$cptfile -T$P -D > $colors;
- m=`echo $P | awk -F"/" 'function max(x,y){return (x>y)?x:y};function abs(x){return (0<x)?x:-x}{print max(abs($1),$2)}'`
- TRIANGLES=`grdinfo -C $U3 | awk -v m=$m 'function max(x,y){return (x>y)?x:y};function abs(x){return (0<x)?x:-x}{if (m<max(abs($6),abs($7))){print "-E"}}'`
+ m=`echo $P | awk -F"/" \
+ 'function max(x,y){return (x>y)?x:y};function abs(x){return (0<x)?x:-x}{print max(abs($1),$2)}'`
fi
my_gmt $INDEX
@@ -374,6 +350,6 @@ while [ "$#" != "0" -o "$Eset" == "1" ];do
fi
fi
- # prevent more empty plots (cancel -E option)
- Eset=""
+ # prevent more empty plots (cancel -O option)
+ Oset=""
done
More information about the CIG-COMMITS
mailing list