1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268
|
#!/bin/csh -f
### ascii2ps -- Generate postscript from 3 column ascii file.
### Example script to generate grd file from phi lambda height
### output of Doris InSAR software.
### Ramon Hanssen, Bert Kampes, 07-Dec-2000.
### a good idea would be to use gmtconvert instead of the ascii input.
### $Id: ascii2ps,v 1.6 2004/09/15 10:49:03 kampes Exp $
### MCC - some GMT v4 compatiblity
################################################################
set PRG = `basename "$0"`
set VER = "v1.1 Doris software"
#set AUT = "bert kampes, (c)1999-2003"
set AUT = "TUDELFT, (c)1999-2009"
echo "$PRG $VER, $AUT"
echo " "
### Handle wrong input
set WRONG = "0"
### Input files: 3 column ascii lon/lat/hei triplets.
set XYZFILE = "$1"
if ( $#argv < 1 ) set WRONG = "1"
if ( $WRONG == "1" ) then
cat << __EOFHD
USAGE:
$PRG 3col_asciifile
Create a postscript file of the data in the 3 column ASCII input file
containing "lon lat hei" triplets (output of lonlathei2ascii
program). GMT is used to interpolate these irregularly distributed
data to a grid.
This script is intended to be transparent and understandable.
Please adapt this script to your own needs. (lambda=longitude).
For example, here one could easily set the correct approximate
height of the terrain.
The output is a postscript file named "temp.ps".
The steps foreseen to be taken are:
[before running this, create a three column ascii file with lon lat hei
using philamh2ascii program.
AND optionally convert this 3 column file to UTM coordinates
using ascii2ascii_UTM
]
2) Call GMT program xyz2grd to perform interpolation of this data.
the number of points in the output grid is 300.
3) call GMT program cpt to create a color table
4) call GMT program to create a postscript file
5) call ghostview to view output.
This script can serve as an example. If you want to have another resolution
output grid, for example exactly 100x100 meters, then you will have
to change this script.
With the public domain software PROJ.4, you can easily convert the
WGS coordinates (Doris output) to X,Y values, and then use this script.
See the script ascii2ascii_UTM for an example how to do this for WGS
to UTM coordinates conversion.
DEPENDENCIES:
GMT tools (gmtset, minmax, xyz2grd, grd2cpt, psbasemap
psscale, grdgradient, psmask, grdimage, pscoast)
standard UNIX (echo, csh, rm, etc.)
ghostview
EXAMPLE:
$PRG lonlathei.dat
$PRG xyz_UTM.dat
SEE ALSO:
cpxfiddle, http://gmt.soest.hawaii.edu/
__EOFHD
exit 1
endif
### Program variables.
set CMAP = sealand
set CPTFILE = temp.cpt
set GRDFILE = temp.grd
set PSFILE = temp.eps
#use paper=a4+ for eps, a4 for ps
rm -f $PSFILE
rm -f $CPTFILE $GRDFILE
#---------------------------
# GMT defaults override
#---------------------------
gmtset MEASURE_UNIT cm \
PAPER_MEDIA a4+ \
PAGE_ORIENTATION portrait \
ANOT_FONT Helvetica \
ANOT_FONT_SIZE 10p \
ANOT_OFFSET 0.2c \
BASEMAP_AXES WeSn \
LABEL_FONT Helvetica \
LABEL_FONT_SIZE 10 \
UNIX_TIME_POS -2c/-2c
#---------------------------
#### program philamh2ascii replaces code below
#### if file temp.xyz exists already, maybe you do not want to wait again?
#### but for now, always recompute it.
#if ( "1" == "2" ) then
#### Convert binary files to ascii row; trick cpxfiddle with -w1.
#echo " Convert LONGITUDES in binary real4 file" \"${LON}\" "to ascii table"
# #cpxfiddle -w1 -qnormal -fcr4 $LON | awk '{printf "%3.4f\n%3.4f\n",$1,$2}' > qlon
# # using new r4 option to cpxfiddle
# cpxfiddle -w1 -qnormal -fr4 $LON > qlon
#echo " Convert LATITUDES in binary real4 file" \"${LAT}\" "to ascii table"
# cpxfiddle -w1 -qnormal -fr4 $LAT > qlat
#echo " Convert HEIGHTS in binary real4 file" \"${HEI}\" "to ascii table"
# cpxfiddle -w1 -qnormal -fr4 $HEI > qhei
#
#### Create xyz file, delete temp files.
#echo "Paste LAT/LON/HEI to xyzfile" \"${XYZFILE}\"
# paste qlon qlat qhei > $XYZFILE.tmp
#echo "removing -999 from input (unwrap problem, or other problem)"
# grep -v '\-999' $XYZFILE.tmp > $XYZFILE
# rm -f qlon qlat qhei $XYZFILE.tmp
#endif
### Get min/max and increment from data itselves.
### Assume lat/lon is formatted as %7.4f?
### Not robust yet?
echo "finding out minmax in ascii column data"
set minmax = `minmax $XYZFILE`
echo minmax $minmax
set WEST = `echo $minmax | cut -d '<' -f2 | cut -d'/' -f1`
set EAST = `echo $minmax | cut -d '<' -f2 | cut -d'/' -f2 | cut -d'>' -f1`
set SOUTH = `echo $minmax | cut -d '<' -f3 | cut -d'/' -f1`
set NORTH = `echo $minmax | cut -d '<' -f3 | cut -d'/' -f2 | cut -d'>' -f1`
set MINHEI = `echo $minmax | cut -d '<' -f4 | cut -d'/' -f1`
set MAXHEI = `echo $minmax | cut -d '<' -f4 | cut -d'/' -f2 | cut -d'>' -f1`
echo " "
### handle 1e+06 format gracefully?
#set EAST = `echo $EAST | awk '{printf "%.20f", $1}'`
#set WEST = `echo $WEST | awk '{printf "%.20f", $1}'`
#set SOUTH = `echo $SOUTH | awk '{printf "%.20f", $1}'`
#set NORTH = `echo $NORTH | awk '{printf "%.20f", $1}'`
echo "Using the following parameters:"
echo "-------------------------------"
echo "WEST: $WEST"
echo "EAST: $EAST"
echo "SOUTH: $SOUTH"
echo "NORTH: $NORTH"
echo "MINHEI: $MINHEI"
echo "MAXHEI: $MAXHEI"
### try to figure out if we have UTM coordinates
### assume utm if north is big.. or small
set CARTESIAN = `echo "$NORTH $SOUTH" | awk '{if(($1>200)||($2<-200)){print"1"}else{print"0"}}'`
echo "CARTESIAN = $CARTESIAN"
#
# numgrd should be multiple of 100?
# number of points in output grid in 1 direction
set NUMGRD = "200"
#set IX = `echo $WEST $EAST $NUMGRD | awk '{print ($2-$1)/($3-1)}'`
#set IX = `echo $WEST $EAST $NUMGRD | awk '{print ($2-$1)/$3}'`
#set IY = `echo $SOUTH $NORTH $NUMGRD | awk '{print ($2-$1)/$3}'`
### report of Bruno Crippa that awk is not precise enough, use bc:
#%// BK 24-Nov-2002
#set IX = `echo "scale=20; ( $EAST - $WEST ) / $NUMGRD" | bc -l`
#set IY = `echo "scale=20; ( $NORTH - $SOUTH ) / $NUMGRD" | bc -l`
# but why not? (awk handles at least 1.2342e+06)
set IX = `echo $WEST $EAST $NUMGRD | awk '{printf "%.20f", ($2-$1)/$3}'`
set IY = `echo $SOUTH $NORTH $NUMGRD | awk '{printf "%.20f", ($2-$1)/$3}'`
#
#set RANGE = `minmax -I$IX/$IY $XYZFILE`
set RANGE = "-R$WEST/$EAST/$SOUTH/$NORTH"
set IXkm = `echo "scale=3; $IX * 40000 / 360" | bc -l`
set IYkm = `echo "scale=3; $IY * 40000 / 360" | bc -l`
echo "RANGE: " $RANGE
if ( "$CARTESIAN" == "1" ) then
echo "Sample interval X: $IX degrees, or $IXkm km at equator"
echo "Sample interval Y: $IY degrees, or $IYkm km at equator"
else
echo "Sample interval X: $IX m"
echo "Sample interval Y: $IY m"
endif
### OR, use minmax to get range, much better...
#set IX = 200
#set IY = 200
#set RANGE = `minmax -I$IX $XYZFILE`
#echo "IX = $IX [m]"
#echo "IY = $IY [m]"
#echo "RANGE = $RANGE"
### Create grd file with xyz2grd (better use surface).
echo "do xyz2grd"
xyz2grd $XYZFILE -G$GRDFILE -I$IX/$IY $RANGE
#awk '{print $1, $2, $3+1000.0}' $XYZFILE | surface -G$GRDFILE -I$IX/$IY $RANGE -T0.75
#surface $XYZFILE -G$GRDFILE -I$IX/$IY $RANGE -T0.75
#echo "for unwrapped ifg, converting to cm slant range displacement"
#grdmath $GRDFILE 1.4 MUL PI DIV = temp.grd2; mv temp.grd2 $GRDFILE
#echo "setting min_z --> 0"
#echo "you may not want this"
# grdmath $GRDFILE $MINHEI SUB = temp.grd2; mv temp.grd2 $GRDFILE
#set MINHEI = 0.0
#set MAXHEI = `echo $MINHEI $MAXHEI | awk '{printf "%f", $2 - $1}'`
### Create colormap from data with grd2cpt.
#set CPTFILE = /home/fmr/d4/hanssen/gmt_scripts/colortables/minpipi.cpt
#set CPTFILE = /home/fmr/d4/hanssen/gmt_scripts/colortables/mass.cpt
echo "creating colormapping..."
set INCHEI = `echo $MINHEI $MAXHEI 255 | awk '{print($2-$1)/$3}'`
echo "grd2cpt $GRDFILE -C$CMAP -Z -S$MINHEI/$MAXHEI/$INCHEI"
#grd2cpt $GRDFILE -C$CMAP -Z -S$MINHEI/$MAXHEI/$INCHEI >! $CPTFILE
grd2cpt $GRDFILE -C$CMAP -Z >! $CPTFILE
### Actually start GMT plotting.
# proj better use correct scaling for x/y
set proj = "-JX15"
set ticks = "-Bf0.1a1.0/f0.1a1.0"
if ( "$CARTESIAN" == "1" ) set ticks = "-Bf5000.0a40000.0/f5000.0a40000.0"
echo "proj = $proj"
echo "ticks = $ticks"
echo "psbasemap $proj $RANGE $ticks -U -K "
psbasemap $proj $RANGE $ticks -U -K >! $PSFILE
echo "psscale -I1 -D12.5/-1.2/5/0.4h -C$CPTFILE -B1000:a:/:m: -O -K "
psscale -I1 -D12.5/-1.2/5/0.4h -C$CPTFILE -B1000:a:/:m: -O -K >> $PSFILE
echo "grdgradient"
grdgradient $GRDFILE -A0 -Nt -G$GRDFILE.gradient
echo "skipping psmask"
#echo "psmask"
#psmask $XYZFILE -I$IX $proj $RANGE -K -O -S100 >> $PSFILE
echo "grdimage"
grdimage $GRDFILE $proj -I$GRDFILE.gradient -C$CPTFILE -R -B -K -O >> $PSFILE
echo "pscoast (last call)"
pscoast $proj -R -B -Dh -W1 -S0/0/255 -O >> $PSFILE
############################################################################
### View result
#if ( $PSFILE > 1kb ) ...
#set GV = ghostview
set GV = gv
echo "using gv to view file: $PSFILE"
#echo "using ghostview to view file: $PSFILE"
$GV $PSFILE
### EOF.
|