Download
1 / 24
240 likes | 399 Views
grdvolume. example18 by MK.Chen Tw-vol by heaven. Auther: 戴富隆. #! /bin/csh ###### File:Geoedu/home/heaven/gmtexp/w05/job18.gmt ###### Author: FL Tai, 11 Nov 2004 ###### GMT EXAMPLE 18 ###### $Id: job18.csh,v 1.1.1.1 2000/12/28 01:23:45 gmt Exp $
E N D
grdvolume example18 by MK.Chen Tw-vol by heaven Auther: 戴富隆
#! /bin/csh ###### File:Geoedu/home/heaven/gmtexp/w05/job18.gmt ###### Author: FL Tai, 11 Nov 2004 ###### GMT EXAMPLE 18 ###### $Id: job18.csh,v 1.1.1.1 2000/12/28 01:23:45 gmt Exp $ ###### Purpose: Illustrates volumes of grids inside contours and spatial ###### GMT progs: gmtset, gmtselect, grdclip, grdcontour, grdgradient, grdimage, ###### grdmath, grdvolume, makecpt, pscoast, psscale, pstext, psxy set inxyz1 = ship.xyz set outgrd1 = AK_gulf_grav.grd set outgrd2 = AK_gulf_grav_i.grd set outgrd3 = mask.grd set outgrd4 = mask-G.grd set outgrd5 = tmp.grd set outps1 = example_18.ps
set output1 = pratt.d set output2 = center.awk set output3 = centers.d set outcpt1 = grav.cpt set range1 = -149/-135/52.5/58 set proj1 = M5.5i set frame1 = 2f1 set frame2 = 2f1WSEn ########################################################################################
# Get Sandwell/Smith gravity for the region #img2latlongrd world_grav.img.7.2 -R$range1 -G$outgrd1 -T1 -V # Use spherical projection since SS data define on sphere #gmtset ELLIPSOID Sphere # Define location of Pratt seamount echo "-142.65 56.25" > $output1 # First generate gravity image w/ shading, label Pratt, and draw a circle # of radius = 200 km centered on Pratt. makecpt -Crainbow -T-60/60/10 -Z -V > $outcpt1 grdgradient $outgrd1 -Nt1 -A45 -G$outgrd2 -V
grdimage $outgrd1 -I$outgrd2 -J$proj1 -C$outcpt1 -B$frame1 \ -P -X1.5i -Y5.85i -K -V > $outps1 pscoast-R$range1 -JM -Di -G160 -W0.25p -O -K -V >> $outps1 psscale -D2.75i/-0.4i/4i/0.15ih -C$outcpt1 -B20f10/:mGal: -O -K -V >> $outps1 $AWK '{print $1, $2, 12, 0, 1, "LB", "Pratt"}' $output1 | \ pstext -R -JM -O -K -D0.1i/0.1i -V >> $outps1 $AWK '{print $1, $2, 0, 200, 200}' $output1 | \ psxy -R -JM -SE -W0.25p -O -K -V >> $outps1 # Then draw 10 mGal contours and overlay 50 mGal contour in green 寫pratt於圖上 畫半徑200km之圓於圖上
grdcontour $outgrd1 -JM -C20 -B$frame2 -U/-1.25i/-0.75i/"Example 18 in Cookbook" \ -Y-4.85i -O -K -V >> $outps1 grdcontour $outgrd1 -JM -C10 -L49/51 -Dsm -Wc0.75p/0/255/0 -O -K -V >> $outps1 pscoast -R -JM -Di -G160 -W0.25p -O -K -V >> $outps1 awk '{print $1, $2, 0, 200, 200}' $output1 | \ psxy -R -JM -SE -W0.25p -O -K -V >> $outps1 \rm -f sm_*[0-9].xyz # Only consider closed contours # Now determine centers of each enclosed seamount > 50 mGal but only plot # the ones within 200 km of Pratt seamount. dumpfile 每個50封閉曲線各自產生sm*檔 封閉曲線上下極限 force. 略過不存在檔案
# First make a simple $AWK script that returns the average position # of a file with coordinates x, y (remember to escape the $ sign) cat << EOF > $output2 BEGIN { x = 0 y = 0 n = 0 } { x += \$1 y += \$2 n++ } END { print x/n, y/n } EOF
# Now determine mean location of each closed contour and # add it to the file centers.d \rm -f centers.d foreach file (sm_*.xyz) awk -f $output2 $file >> $output3 end # Only plot the ones within 200 km gmtselect -R -JM -C200/$output1 $output3 -V > $$ psxy $$ -R -JM -SC0.04i -G255/0/0 -W0.25p -O -K -V >> $outps1 psxy -R -JM -ST0.1i -G255/255/0 -W0.25p $output1 -O -K -V >> $outps1 用output2 迴圈平均file(sm_*.xyz) seamount 位置 Progfile 搜索檔案 subset Radius/center
# Then report the volume and area of these seamounts only # by masking out data outside the 200 km-radius circle # and then evaluate area/volume for the 50 mGal contour grdmath -R -I2m -F -142.65 56.25 GDIST = $outgrd3 grdclip $outgrd3 -Sa200/NaN -Sb200/1 -V -G$outgrd4 grdmath $outgrd1 $outgrd4 -V MUL = $outgrd5 set info = `grdvolume $outgrd5 -C50 -Sk` Outgrd1*outgrd4 終端機執行:4 column 50 14157.7 382970 27.0502 Z area vol vol/area
psxy -R -JM -A -L -W0.75p -G255 -O -K -V << EOF >> $outps1 -148.5 52.75 -140.5 52.75 -140.5 53.75 -148.5 53.75 EOF pstext -R -JM -O -V << EOF >> $outps1 -148 53.08 14 0 1 LM Areas: $info[2] km@+2@+ -148 53.42 14 0 1 LM Volumes: $info[3] mGal\264km@+2@+ EOF ######################################################################################## #\rm -f $$ grav.cpt sm_*.xyz *_i.grd tmp.grd mask.grd pratt.d center* .gmt* #imagetool $outps1 & gv $outps1 &
Sm_50_6i.xyz Set value50 • 223.483 57.7991 50 • 223.517 57.7976 50 • 223.55 57.7878 50 • 223.558 57.7833 50 • 223.583 57.7665 50 • 223.6 57.75 50 • 223.61 57.7167 50 • 223.604 57.6833 50 • 223.583 57.661 50 • 223.553 57.65 50 • 223.55 57.6496 50 • 223.548 57.65 50 • 223.517 57.6612 50 • 223.49 57.6833 50 • 223.483 57.6908 50 • 223.467 57.7167 50 • 223.455 57.75 50 • 223.464 57.7833 50 • 223.483 57.7991 50 colse
Output2 center.awk BEGIN { x = 0 y = 0 n = 0 } { x += $1 y += $2 n++ } END { print x/n, y/n }
Output3 center.d • 213.641 56.4759 • 216.948 56.4144 • 214.169 56.3746 • 214.735 56.3001 • 217.35 56.2427 • 215.663 56.044 • 216.786 55.9797 • 218.128 55.8274 • 217.197 55.4877 • 219.637 55.1186 • 216.906 54.9629 • 223.029 54.5248 • 212.767 54.4857 • 211.494 54.103 • 211.646 54.4714 • 222.596 54.0813 • 212.214 54.1027 • 212.11 53.7308 • 223.394 53.6759 • 215.67 53.5724 • 223.901 53.4917 • 224.364 53.3184 • 215.913 53.0323 • 223.531 57.7252 • 223.836 57.3779 • 223.916 57.1381 • 217.362 57.0014
grdcontour contouring of 2-D gridded data sets grdcontourgrdfile-Ccont_int-Jparameters [ -A[- ][anot_int][ffont_size][aangle] [/r/g/b][o]] ] [ -Btickinfo ][ -Ddumpfile ] [ -Eazimuth/elevation ] [ -Ggap/width ] [ -K ] [ -Llow/high ] [ -M[flag] ] [ -N[[-]unit] ] [ -O ] [ -P ] [ -Qcut ] [ -Rwest/east/south/north[r] ] [ -Ssmoothfactor ] [ -T[+|-][gap/length][:LH] ] [ -U[/dx/dy/][label] ] [ -V ] [ -W[+][type]pen ] [ -Xx-shift ] [ -Yy-shift ] [ -Z[factor[/shift]][p] ] [ -ccopies ] [ -bo[s] ] -D Dump the (x,y,z) coordinates of each contour to separate files, one for each contour segment. The files will be named dumpfile_cont_segment[_i].xyz (or .b is -b is selected), where cont is the contour value and segment is a running segment number for each contour interval (for closed contours we append _i.) However, when -M is used in conjunction with -D a single mul- tisegment file is created instead.
gmtselect • Select data subsets based on multiple spatial criteria • gmtselect [ infiles ] [ -Amin_area[/min_level/max_level] ] [ -C[f]dist/ptfile ] [ -Dresolution ] [ -Fpolygonfile ] [ -H[nrec] ] [ -I[cflrs] ] [ -Jparameters ] [ -Ldist/linefile ] [ -M[flag] ] [ -Nmaskvalues[o] ] [ -Rwest/east/south/north[r] ] [ -V ] [ -: ] [ -bi[s][n] ] [ -bo[s] ] • -C Pass all records whose location is within dist km of any of the points in ptfile. If dist is zero then the 3rd column of ptfile must have each point's individual radius of influence. Prepend f to indicate you want approximate flat earth distance calculations (faster) than exact great circle calculations (slower). • -S Convert degrees to meters, append k for km [Default is Cartesian].
grdmath • Reverse Polish Notation calculator for grd files • grdmath [ -Ixinc[m|c][/yinc[m|c]] -Rwest/east/south/north -V] operand [ operand ] OPERATOR [ operand ] OPERATOR ... = outgrdfile • -I x_inc [and optionally y_inc] is the grid spacing. Append m to indicate minutes or c to indicate seconds. • GDIST 2 Great distance (in degrees) between grid nodes and stack lon,lat.
grdvolume Calculating volume under a surface within a con- tour grdvolumegrdfile [ -Ccval or -Clow/high/delta ] [ -Lbase ] [ -Rwest/east/south/north[r] ] [ -S[k] ] [ -T ] [ -V[l] ] [ -Zfact[/delta] ] -C find area and volume inside the cval contour. Alterna- tively, search using all contours from low to high in steps of delta. [Default returns entire area and volume of grid].
