grd2xyz 解 etopo2 (netcdf) 然後用gawk挖空其中一塊…

自從GMT支援netcdf 後，有些三維的nc檔可以透過grd2xyz解…還算方便。
下面的程式主要是我有二個水深檔，一個是台灣附近500m解析度，一個是etopo2。
我想取出東經114.5-150.5，北緯14.5-55.5的水深，有500m解析度的水深就用他，沒有的話就用etopo2。然後合起來的水深內插成5分一個點，並轉成xyz檔。程式如下：

R1="-R114.5/150.5/14.5/55.5"
slonl=117
slonr=125
slatt=27
slatb=18
FN500m=taidpv626_500m.xyz
OFN=decetop.dat
DELTX="-I5m/5m -S2m -N1"
grd2xyz ETOPO2v2g_f4.nc -fg $R1 | gawk '! ($1 > L && $1 < R && $2 > B && $2 < T) {printf"%8.4f %8.4f %7.1f\n", $1, $2, -1*$3}' L=$slonl R=$slonr T=$slatt B=$slatb > $OFN
gawk '{print $1, $2}' decetop.dat | psxy $R1 -JM5i -S+0.01 -W0.001c -Ba10f5SWne -P -K > decetop.ps
pscoast -R -JM -B -G104 -Df -O -W0.01,104 >> decetop.ps
cat $FN500m >> $OFN
makecpt -Cseis -T-8000/8000/500 -Z > colors.cpt
nearneighbor $OFN $R1 $DELTX -Gbath.grd 1>/dev/null 2>&1
grdview bath.grd $R1 -JM5i -Qi100 -V -Sc -Ccolors.cpt -P -K > bath.ps
pscoast -R -JM -B -G104 -Df -O -W0.01,104 >> bath.ps
grd2xyz bath.grd -fg $R1 | gawk  '{printf"%8.4f %8.4f %7.1f\n", $1, $2, $3}'> resample_$OFN

ps2epsi decetop.ps decetop.eps
eps2jpg -d 600 -f decetop.eps
ps2epsi bath.ps bath.eps
eps2jpg -d 600 -f bath.eps
rm -f decetop.ps decetop.eps bath.ps bath.eps