GMT学习总结
1、了解投影法在使用GMT前，必需先瞭解所使用的投影法。
地球的幾何模型可以被投影到多種表面上，常用的有
(1)方位投影(2)圓錐投影(3)圓柱投影。
一般畫地形圖通常使用「橫麥卡托投影法」(Mercator)，屬於圓柱投影的一種
，將全球按經度分為60個6度分帶，在GMT裡表示的方法為-Jm或-JM。而GMT
內定的GPS座標系統（ELLIPSOID）則使用WGS84座標系統。
-Jm6i代表在紙上每一度畫6英吋；-JM6i代表一個圖，不論幾度，畫出來的圖，
總長6英吋，如果要改成公分，請加上參數c，如-Jm6c。
2、開始繪圖之前，先瞭解GMT預設的參數值（gmtdefaults）
GMT有許多預設的設定，要看GMT預設參數，在unix環境下打 gmtdefaults -D。
若要看使用者現正使用的參數，打-L。
預在job內改參數，gmtset + 要改的參數，例如：
gmtset ANOT_FONT_SIZE 12 （字型大小改成12）
gmtset PAPER_MEDIA A4 (紙張大小設成A4，若要畫A0或A1的大圖，則改為A0或A1)
gmtset BASEMAP_TYPE FANCY (畫出的地形圖邊框為一黑一白間格，此為GMT預設值)
3、瞭解水深資料格式
本實驗室所使用之水深資料，大部份來自國家海洋科學研究中心海洋資料庫。
資料來源有：
(1) 海研一號、二號、三號研究船多年來所收集之船測水深資料
(2) MEST、ETOPO5、SCSBDv1、GTOPO30、TaiDBMv6等網格化資料
(3) 其他人工數畫海圖
綜合以上資料整編而來，對水深資料有興趣者，請參考國家海洋科學
中心_水深資料庫網頁，http://www.ncor.ntu.edu.tw/ODBS/depth/ek500.html。
常用的台灣週圍海域數值地形模型(DEM)為netCDF數值格式化檔案
（taidp.v624.grdint_500m.grd），其解析度為每500公尺一點，包含
的範圍自北緯18度到27度、東經117度到125度。欲查詢grd檔包含的
範圍，請打grdinfo *.grd。
4、常用的GMT參數
-B
畫圖的外框，常用的參數有：a(annotation)、f(frame)、g(grid)，例如：-Ba10f5g2即每十度加註一個標註、每五度畫一個TICK、每2度畫一條線，另外，GMT的預設值是WESN，即畫框的東南西北邊都會標示註解和TICK，若只要畫東邊和北邊，則寫成-Ba10f5g2wEsN(要畫出來的那邊用大寫)。
-J
選擇投影法，如-Jm6i
-R
Range，依投影法的不同，指定畫圖的範圍，如-R20/23/119/120，就是畫出北緯20度到23度、東經119度到120度
-U
-U ]上標籤的位置跟標籤文字
-K及-O
這兩個參數是剛開始使用GMT時，很容易搞混的地方，-K與-O控制多層圖檔的疊加，-K代表後來的圖層可以疊加在這個圖層上，-O是代表Overlay。在script的第一個程式裡，一定要加-K，接下來的第二、第三…到第N-1個程式裡，若是要疊加在前一個圖上，全部都要加-K及-O，在第N個程式時（即最後一個），加-O即可
-P
GMT預設的紙張是橫放，加了-P即變成直放
-X、-Y
設定從X軸及Y軸的何處當成起始點。GMT的預設值是以紙張的左下角(0,0)為起始點，若改成-X5、-Y6即從(5,6)開始當成畫圖的起始點
-:
GMT讀入的檔案，通常第一欄是經度，第二欄是緯度，加了-:這是把兩欄的位置換過來，故若一個檔是緯度在前，經度在後，即須要用-:
5、如何畫3D地形 grdview
範例：3d_tech.bash
#!/bin/bash
# GMT V3.4
scale=6i
zscale=0.0002
cptf=taidp2.cpt
proj=M
orient=-P
gmtset ANOT_FONT_SIZE 12
gmtset PAPER_MEDIA A4
fram=a1f30m
inputgrdf=taidp.v624.grdint_500m.grd
range=121/123/23/25:30
outputf=3d_180_90_teach.ps
viewport=180/90
#### topo_3d
grdgradient $inputgrdf -A30 -Ggrad.grd -V
grdhisteq grad.grd -Ggradhisteq.grd -V -N
grdmath gradhisteq.grd 0.5 x = intensity.grd
grdview inputgrdf−Rrange/-10000/4000 -Bfram −Jprojscale−Jzzscale -Qi100 -V -K -Ccptf orient -Eviewport−Iintensity.grd −N−10000/255/255/255>outputf
scale=6i
設定比例尺大小
zscale=0.0002
設定z方向的比例尺大小
cptf=taidp2.cpt
設定使用到的cpt檔
proj=M
設定投影法種類，M為橫麥卡托投影法
orient=-P
設定紙張的方向，-P為直放
gmtset ANOT_FONT_SIZE 12
設定邊框註記字體的大小為12
gmtset PAPER_MEDIA A4
設定紙為A4 size
fram=a1f30m
設定邊框每1度標一數字，每30分畫一個tickmark
inputgrdf=taidp.v624.grdint_500m.grd
輸入的水深檔
range=121/123/23/25:30
設定圖的範圍
outputf=3d_180_90_teach.ps
設定寫出的檔案名稱
viewport=180/90
方位角180度(azimuth)，仰角90度(elevation)，看起來像水平2D圖
6、如何畫等深線 grdcontour
#!/bin/bash
# GMT V3.4
#
###### general setting
scale=6i
proj=M
orient=-P
gmtset DEGREE_FORMAT 3
gmtset ANOT_FONT_SIZE 12
gmtset PAPER_MEDIA A4
gmtset BASEMAP_TYPE FANCY
inputf=taidp.v624.grdint_500m.grd
range=119/120:30/21:30/23
outputf=contour_tech.ps
###### contour
grdcontour inputf−Rrange -Jprojscale\
orient−Ba1f30mEwnS−L−10000/0 −C200−Wa2/0/0/250−Wc1/0−A1000f5 −K>outputf
7、如何產生彩色CPT檔(1) grd2cpt
方法一：grd2cpt (根據grd檔的水深值，來產生cpt檔)
grd2cpt taidp.v624.grdint_500m.grd -Crainbow -Z -S-6000/4000/500 >teach_01.cpt
grd2cpt taidp.v624.grdint_500m.grd -Crainbow -S-6000/4000/500 >teach_02.cpt
範例：image_tech.bash
#!/bin/bash
# GMT V3.4
scale=4i
proj=M
orient=-P
inputf=taidp.v624.grdint_500m.grd
range=121/124/22:30/25
cptf1=teach_01.cpt
cptf2=teach_02.cpt
outputf=image_tech.ps
##### teach_01.cpt
grdimage inputf−Rrange -V -Jprojscale -Ba1f30mEwnS:"teach_01.cpt":N orient−K −Ccptf1 >outputfpscoast−Rrange orient−Jprojscale−Df −W4/0−V−K−O>>outputf
psscale -Ccptf1−D2/−0.4/4/0.2h −K−Oorient -V >> outputfpstext−Jprojscale−Vorient -O -K -R -V -G0 <<END>> outputf 123 24:50 12 0 1 1 teach_01.cpt(-Z) END ##### teach_02.cpt grdimageinputf -Rrange−V−Jprojscale−Y5.3−Ba1f30mEwnSorient -K -O\
-Ccptf2>>outputf
pscoast -Rrangeorient -Jprojscale -Df \
-W4/0 -V -K -O >>outputfpsscale−Ccptf2 -D2/-0.4/4/0.2h \
-K -O orient−V>>outputf
pstext -Jprojscale -V orient−O−K−R−V−G0<<END>>outputf
123 24:50 12 0 1 1 teach_02.cpt
END
8、如何產生彩色CPT檔(2) makecpt
方法二：makecpt (直接產生cpt檔)
makecpt -Crainbow -T-7000/4000/500>makecpt_01.cpt
9、如何加站位 psxy
範例：psxy_tech.bash
#!/bin/bash
# GMT V3.4
scale=4i
zscale=0.0002
cptf=taidp2.cpt
inputf=taidp.v624.grdint_500m.grd
proj=M
orient=-P
gmtset ANOT_FONT_SIZE 10
fram=a30mf15mwENs
range=119/120:45/20:55/22:30
outputf=psxy_tech.ps
####################
# 畫等深線
grdcontour inputf−Rrange \
-Jprojscale orient−Bfram -L-10000/0\
-C500 -Wa2/0 -Wc1/50 -A1000f5 -K >$outputf
####################
# 加站位(紫色三角形)
psxy 0421A.dat -Jprojscale \
-K -O orient−St0.1 −G250/0/250−Rrange -V>>$outputf
#####################
# 畫海岸線
pscoast -Rrangeorient -Jprojscale -Df \
-G200/200/200 -W4/2 -V -O -K>>$outputf
10、如何加字--pstext
範例：pstext_tech.bash
#!/bin/bash
# GMT V3.4
scale=5i
zscale=0.0002
cptf=taidp2.cpt
inputf=taidp.v624.grdint_500m.grd
proj=M
orient=-P
gmtset ANOT_FONT_SIZE 10
gmtset DEGREE_FORMAT 2
fram=a30mf15mwENs
range=119/120:45/20:55/22:30
outputf=pstext_tech.ps
####################
# 加站位點(橘色三角)
psxy 0421A.dat -Jprojscale \
-K orient−St0.1−G255/178/0−Rrange -V>$outputf
####################
# 加站名編號(綠色數字)
pstext 0421A.dat -R orient−V−Jprojscale −K−O−G0/150/0>>outputf
##########
#加註解文字
pstext -Jprojscale -V orient−O−K−R−V−G0 <<END>>outputf
119:15 21:10 16 0 2 1 Demo. Well Site
119:15 21:20 16 0 1 1 Southern Taiwan
END
###########
# 加海岸線
pscoast -Rrangeorient -Jprojscale -Bfram−Df −G200/200/200−W4/2−V−O−K>>outputf
11、如何画比例尺 psbasemap
#!/bin/bash
# GMT V3.4
scale=4i
proj=M
orient=-P
gmtset DEGREE_FORMAT 2
gmtset ANOT_FONT_SIZE 12
gmtset PAPER_MEDIA A4
gmtset BASEMAP_TYPE FANCY
inputf=taidp.v624.grdint_500m.grd
range=121/124/22:30/25
cptf=teach_01.cpt
outputf=psbasemap_tech.ps
############
# 画彩色地形
grdimage inputf−Rrange -V -Jprojscale orient−K −Ccptf >$outputf
############
# 画比例尺
psbasemap -Rrange−Jprojscale −Lf122:30/25.3/25.5/100−V−K−O>>outputf
############
# 画海岸线
pscoast -Rrangeorient -Jprojscale \
-Ba1f30mWSne:."psbasemap": -Df \
-W4/0 -V -K -O>>$outputf
############
# 画colorbar
psscale -Ccptf−D2/−0.4/4/0.2h −K−Oorient -V >> $outputf
12、如何畫研究船軌跡圖 (1) timemark
寫在前面：
研究船的航跡資料(ASCII格式), 乃是每個航次研究船出海作業的路線, 大部份都經劉紹勇學長處理過, 成為gmt容易讀取的格式, 依照各航次的不同分別放在不同的目錄下
最簡單的航線圖就是只把測線畫出來, 上面沒有做任何標記, 但是在做震測資料處理時, 常需要畫出A0的大圖來對照, 此時就需要在測線上做標記, 常用的標記有三種：(1)timemark (2)ffid (3) cdp, 每種資料需經處理才能得到gmt可使用的格式
產生timemark步驟
步驟 1
在dtk20a上, 執行timemark程式
步驟 2
name of input file = ＞(輸入檔名)or1_524b.raw.int
步驟 3
tick mark per ? minutes = = ＞ (幾分鐘標一個tickmark？)30
步驟 4
tick mar per ? hours = = ＞ (幾個小時標一個timemark？)2
步驟 5
size of marks <3-20 is reasonable > = = ＞ (寫出的tickmark字型大小為何？)5
步驟 6
最後產生tickmark.txt 及 timemark.txt檔, 最後由GMT畫出
畫出船行測線圖
範例：timemark_tech.bash
#!/bin/bash
# GMT V3.4
scale=5i
proj=M
orient=-P
gmtset DEGREE_FORMAT 2
gmtset ANOT_FONT_SIZE 12
gmtset PAPER_MEDIA A4
gmtset BASEMAP_TYPE FANCY
inputf=taidp.v624.grdint_500m.grd
range=121:30/124/24:30/26
cptf=teach_01.cpt
outputf=timemark_tech.ps
############
# 標 timemark
pstext 524b_timemark.txt -Rrangeorient -V -Jprojscale -K -G0/0/250 > $outputf
############
# 標tickmark
pstext 524b_tickmark.txt -R orient−V−Jprojscale−K−O−G0/0/250>>outputf
############
# 標測線位置
psxy or1_524b.raw.int -Jprojscale -R -V orient−Sc0.01−O−K−G250/0/250>>outputf
############
# 畫海岸線
pscoast -Rrangeorient -Jprojscale -Ba1f30mWSne:."timemark_teach": -Df -W4/0 -V -K -O>>$outputf
13、台灣東部海域深度20M海流圖
#! /bin/csh
###### Author: fingerfish
###### GMT HOMEWORK2
###### Purpose:
######3-D perspective color plot of Africa topography and img
###### GMT progs: grdcontour, grdview, pscoast, pstext, psxyz
###### DATA FROM GMTDEMO
set ingrd1 = africa_geoid.grd
set ingrd2 = africa_topo.grd
set outgrd1 = g_intens.grd
set outgrd2 = t_intens.grd
set outps = africa2.ps
set incpt1 = geoid.cpt
set incpt2 = topo.cpt
set proj1 = M5i
set proj2 = Z1.3i
set proj3 = x1i
set range = -20/60/-40/40
set frame1 = a10f5/a10f5NESW
set frame2 = 2/2/2000
set viewpoint = 210/30
##############################################################
#set GMT defaults
gmtset DEGREE_FORMAT 3
#to make img files transfer to grd files
img2grd ~/geoid.9.2.img -R$range -Gafricatmp.grd -V -T0
#Resample a grd file onto a new grid
grdsample africatmp.grd -Gingrd1−I10m−Rrange
grdsample ~/topo5.grd -Gingrd2−I10m−Rrange
#create color palette table for topogragpy
makecpt -Crainbow -T-6500/6000/200 -Z > incpt1makecpt−Crelief−T−6000/4000/50−Z>incpt2
#create grd file for illuminnation
grdgradient ingrd1−A0−Goutgrd1 -Nt0.75 -M
grdgradient ingrd2−A0−Goutgrd2 -Nt0.75 -M
#Create 2-D perspective grayshaded/colored image of African geoid
grdview ingrd1−Ioutgrd1 -Rrange−Jproj1 -Eviewpoint−Cincpt1 -Qs \
-U/-1.25i/-1i/"Example 4c in Cookbook" -P -X0.6i -Y1.25i -K -V > $outps
#make contour lines
grdcontour ingrd1−R−Jproj1 -Eviewpoint−C500−W0.1−O−K−V>>outps
#coastline and river database
pscoast -R -JM -Bframe1−Eviewpoint -W2 -N1 -O -K -V >> $outps
#Arrow pointing the nothern direction
echo '-35 0 0 0 1.1' | \
psxyz -R -JM -Eviewpoint−SV0.2i/0.5i/0.4i−W1p−G255/0/0−N−O−K−V>>outps
echo '-38 -8 36 0 1 5 N' | \
pstext -R -JM -Eviewpoint−N−O−K−V>>outps
#Create 3-D perspective grayshaded/colored image of African topography
grdview ingrd2−Ioutgrd2 -JM -Jproj2−Cincpt2
-Eviewpoint−Rrange'/-6000/4000' \
-N-6000/200/200/200 -Qi100 -Y2.2i -O -K -V >> $outps
#To plot 3-D basemaps at Z=-6000m
psbasemap -R -JM -Bframe2′:Topo(m):wsZ′−Jproj2
-Eviewpoint−Z−6000−O−K−V>>outps
#pstext area
echo '3.25 5.75 60 0.0 33 2 @#AFRICA@# @#RIDGE' | \
pstext -R0/10/0/10 -Jx1i -O -V >> $outps
##############################################################
#remove extrafiles
\rm -f *tmp.grd .gmtcommands
#ghostview the outputs
gv $outps &
14、FPS(非洲)
#! /bin/csh
###work01 Africa
###date 7/22/05
###name fingerfish
####################
set ingrd1 = africa.grd
set ingrd2 = africa_intens.grd
set outcpt = AFRICA.cpt
set infps = FPS.dat
set outps = AFRICA.ps
set range = -20/60/-40/40
set proj = M7i
set frame = a10f5/a10f5NEWS
#set GMT defaults
gmtset ANOT_FONT_SIZE 10 LABEL_FONT_SIZE 10
#creat color palette table for topography
makecpt -Crelief -T-6000/6000/50 -Z > $outcpt
#create grd file for illuminnation
grdgradient ingrd1−Gingrd2 -A0 -M -Nt0.75b -V
###Start plotting
#color image of the topography
grdimage ingrd1−Iingrd2 -Rrange−Jproj -Bframe−Coutcpt -X1.5 -Y4.5 -P -K -V > $outps
#display color scale
psscale -Coutcpt−D8.9/−1/17.8/0.3h−B1000:Topography:−O−K−V>>outps
#coastline and river database
pscoast -R -JM -Di -I1/3/0/0/255 -N1/2 -A2 -W3/0/0/255 -K -O -V >> $outps
###contour line
grdcontour ingrd1−R−Jproj -O -K -S -C500 >> $outps
###plot harvart CMT
awk '{print 1,2,3,4,5,6,7,8,9,10,1,2}' infps|psmeca−JM−R−Sm0.2−G255/0/0−O−K−V>>outps
###pstest
#pstext -R -JM -G255/255/255 -N -O -K -V << END >> $outps
#345 -30 50 0 0 1 Africa
#END
#pstext -R -JM -G0/0/0 -N -O -V << END >> $outps
#20 50 50 0 0 CM AFRICA
#20 45 50 0 0 CM TOPOGRAPHY
END
###ghostview
gv $outps&
15、非洲3-D地形圖以及
#! /bin/csh
###### Author: fingerfish
###### GMT HOMEWORK2
###### Purpose: 3-D perspective color plot of Africa topography and img
###### GMT progs: grdcontour, grdview, pscoast, pstext, psxyz
###### DATA FROM GMTDEMO
set ingrd1 = africa_geoid.grd
set ingrd2 = africa_topo.grd
set outgrd1 = g_intens.grd
set outgrd2 = t_intens.grd
set outps = africa2.ps
set incpt1 = geoid.cpt
set incpt2 = topo.cpt
set proj1 = M5i
set proj2 = Z1.3i
set proj3 = x1i
set range = -20/60/-40/40
set frame1 = a10f5/a10f5NESW
set frame2 = 2/2/2000
set viewpoint = 210/30
################################################################
#set GMT defaults
gmtset DEGREE_FORMAT 3
#to make img files transfer to grd files
img2grd ~/geoid.9.2.img -Rrange -Gafricatmp.grd -V -T0 #Resample a grd file onto a new grid grdsample africatmp.grd -Gingrd1 -I10m -Rrangegrdsample /topo5.grd−Gingrd2 -I10m -Rrange #create color palette table for topogragpy makecpt -Crainbow -T-6500/6000/200 -Z >incpt1
makecpt -Crelief -T-6000/4000/50 -Z > incpt2 #create grd file for illuminnation grdgradientingrd1 -A0 -Goutgrd1−Nt0.75−Mgrdgradientingrd2 -A0 -Goutgrd2 -Nt0.75 -M #Create 2-D perspective grayshaded/colored image of African geoid grdviewingrd1 -Ioutgrd1−Rrange -Jproj1−Eviewpoint -Cincpt1−Qs −U/−1.25i/−1i/"Example4cinCookbook"−P−X0.6i−Y1.25i−K−V>outps
#make contour lines
grdcontour ingrd1−R−Jproj1 -Eviewpoint−C500−W0.1−O−K−V>>outps
#coastline and river database
pscoast -R -JM -Bframe1−Eviewpoint -W2 -N1 -O -K -V >> outps #Arrow pointing the nothern direction echo '-35 0 0 0 1.1' | \ psxyz -R -JM -Eviewpoint -SV0.2i/0.5i/0.4i -W1p -G255/0/0 -N -O -K -V >> outpsecho′−38−836015N′| pstext−R−JM−Eviewpoint -N -O -K -V >> outps #Create 3-D perspective grayshaded/colored image of African topography grdviewingrd2 -Ioutgrd2−JM−Jproj2 -Cincpt2−Eviewpoint -Rrange′/−6000/4000′ −N−6000/200/200/200−Qi100−Y2.2i−O−K−V>>outps
#To plot 3-D basemaps at Z=-6000m
psbasemap -R -JM -Bframe2′:Topo(m):wsZ′−Jproj2 -Eviewpoint−Z−6000−O−K−V>>outps
#pstext area
echo '3.25 5.75 60 0.0 33 2 @#AFRICA@# @#RIDGE' | \
pstext -R0/10/0/10 -Jx1i -O -V >> outps ########################################################## #remove extrafiles \rm -f *tmp.grd .gmtcommands #ghostview the outputs gvoutps &
16、MODLE/coastline/image/contour/label(非洲)
#! /bin/csh
###work01 Africa
###Auther fingerfish 7/28/2005
####################
set ingrd1 = africa.grd
set ingrd2 = africa_intens.grd
set outcpt = AFRICA.cpt
set outcpt2 = AFRICA2.cpt
set infps = FPS.dat
set outps = AFRICA.ps
set range = -20/60/-40/40
set proj = C20/0/7i
set frame = a10f5g10/a10f5g10NEWS
#set GMT defaults
gmtset ANOT_FONT_SIZE 10 LABEL_FONT_SIZE 10
#creat color palette table for topography
makecpt -Crelief -T-6000/6000/20 -Z > outcptmakecpt−Crainbow−T−6000/6000/2000−Z>outcpt2
#create grd file for illuminnation
grdgradient ingrd1−Gingrd2 -A0 -M -Nt0.75b -V
###Start plotting
#color image of the topography
grdimage ingrd1−Iingrd2 -Rrange−Jproj -Bframe−Coutcpt -X1.5 -Y4.5 -P -K -V > $outps
#display color scale
psscale -Coutcpt−D8.9/−1/17.8/0.3h−B1000:Topography:−O−K−V>>outps
#coastline and river database
pscoast -R -Jproj−Di−I1/3/0/0/255−N1/2−A2−W3/0/0/255−K−O−V>>outps
###contour line
grdcontour ingrd1−Coutcpt2 -A2f6 -R -Jproj−O−K−S−C1000>>outps
###plot harvart CMT
awk '{print 1,2,3,4,5,6,7,8,9,10,1,2}' infps|psmeca−Jproj -R -Sm0.2 -G255/0/0 -O -K -V >> $outps
###pstest
pstext -R -Ja -G255/255/255 -N -O -V << END >> $outps
10 10 50 0 0 CB AFRICA
END
###ghostview
gv $outps&
17、Fault-plane solutions,coastline
#! /bin/csh
###### Author:benjones
###### http://www.seismology.harvard.edu
###########################################################
#set gmt defaults
gmtset ANOT_FONT_SIZE 10 LABEL_FONT_SIZE 10
#create color palette table for topography
makecpt -Crainbow -T-8100/4000/50 -Z -V > t12.cpt
#create grd file for illumination
grdgradient t12.grd -Gt12out.grd -A0 -M -Nt0.75 -V
#color image of the topography
grdimage t12.grd -It12out.grd -R-80/10/0/60 -JM14 -Ba10f10/a10f10NEWS -Ct12.cpt -X6 -Y5 -K -V >t12.ps
#add length scales
psbasemap -R -JM14 -B10g10 -K -O -V -Lx4/6/-40/600 >>t12.ps
#display color scale
psscale -Ct12.cpt -B2000 -D7/-1/13/0.1h -K -O -V >>t12.ps
#draw all rivers and national borders and coastline
pscoast -R -JM -Dh -Ia/2/0/0/255 -A2 -W2/150/150/150 -K -Na/3/255/0/255 -O -V >>t12.ps
#plot fault-plane solutions
awk '{print 1,2,3,4,5,6,7,8,9,10,1,2}' fps.dat | psmeca -JM -R -Sm0.2 -G255/0/0 -O -K -V >> t12.ps
#add symbol
psxy -JM -R -Sc0.1 -G255/0/0 -N -O -K -V << END >>t12.ps
325 30
END
#add text
pstext -R -JM -G0/0/0 -N -O -V << END >>t12.ps
325 65 20 0 10 CM Northern Atlantic Ocean Earthquakes
END
18、volume and spatial selections
#! /bin/csh
#grdsample grav.11.2.grd -Gna_grav.grd -R-35/-24/18/25 -I0.5m
###### Author: benjones yang, 2005
###### GMT EXAMPLE 18
###### 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 = na_grav.grd
set outgrd2 = na_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 = -35/-24/18/25
set proj1 = M5.5i
set frame1 = 2f1
set frame2 = 2f1WSEn
########################################################################################
# Get Sandwell/Smith gravity for the region
#img2latlongrd world_grav.img.7.2 -Rrange1−Goutgrd1 -T1 -V
# Use spherical projection since SS data define on sphere
gmtset ELLIPSOID Sphere
# Define location of Pratt seamount
echo "-31 21" > output1 # First generate gravity image w/ shading, label Pratt, and draw a circle # of radius = 200 km centered on Pratt. makecpt -Crainbow -T-400/300/10 -Z -V >outcpt1
grdgradient outgrd1−Nt1−A45−Goutgrd2 -V
grdimage outgrd1−Ioutgrd2 -Jproj1−Coutcpt1 -Bframe1 −P−X1.5i−Y5.85i−K−V>outps1
pscoast -Rrange1−JM−Di−G160−W0.25p−O−K−V>>outps1
psscale -D2.75i/-0.4i/4i/0.15ih -Coutcpt1−B100f10/:mGal:−O−K−V>>outps1
awk '{print 1,2, 12, 0, 1, "LB", "hav"}' 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
grdcontour outgrd1−JM−R−35/−24/18/25−C200−Bframe2 -U/-1.25i/-0.75i/"Example 18 in Cookbook" \
-Y-4.85i -O -K -V >> outps1grdcontouroutgrd1 -JM -R-35/-24/18/25 -C10 -L299/301 -Dsm -Wc0.75p/0/255/0 -O -K -V >> outps1pscoast−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_*.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.
# 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 -foutput2 file>>output3
end
# Only plot the ones within 200 km
gmtselect -R -JM -C200/output1output3 -V >
psxy
-R -JM -SC0.04i -G255/0/0 -W0.25p -O -K -V >> outps1psxy−R−JM−ST0.1i−G255/255/0−W0.25poutput1 -O -K -V >> outps1 # 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 -I0.5m -F -329 21 GDIST =outgrd3
grdclip outgrd3−Sa200/NaN−Sb200/1−V−Goutgrd4
grdmath outgrd1outgrd4 -V MUL = outgrd5setinfo=`grdvolumeoutgrd5 -C200 -Sk`
psxy -R -JM -A -L -W0.75p -G255 -O -K -V << EOF >> outps1329.5 18.2335.918.2335.919.2329.519.2EOFpstext−R−JM−O−V<<EOF>>outps1
329.7 18.45 14 0 1 LM Areas: infokm@+2@+329.718.951401LMVolumes:info mGal\264km@+2@+
EOF
########################################################################################
\rm -f $$ grav.cpt sm_*.xyz *_i.grd tmp.grd mask.grd pratt.d center* .gmt*
#imagetool $outps1 &
gv $outps1 &
19、Colored image of topography, scale, coastline and map frame
#! /bin/csh
###### homework1.gmt
###### Author : didida 10/23/2006
###### purpose : Clored image of topography,scale,coastlin eand map frame
set worldgrd : TOPO5.GRD
set ingrd1 : China.grd
set ingrd2 : China_intens.grd
set outcpt : China.cpt
set outps : China.ps
set range : 72/136/17.5/60
set proj : M7i
set frame : a10f5/a5f5neWS
###########################################################
#set GMT defaults
gmtset ANOT_FONT_SIZE 10 LABEL_FONT_SIZE 15
#Creat color palette table for topography
makecpt -Chaxby -T-7600/8800/100 -Z > outcpt #grdsample grdsampleworldgrd -Gingrd1−Rrange -V -I10m
#Create grd file for illumination
grdgradient ingrd1−Gingrd2 -A0 -M -Nt0.75b -V
###Start plotting
#Color image of the topography
grdimage ingrd1−Iingrd2 -Rrange−Jproj -Bframe−Coutcpt -X1 -Y2 -K -V > outps #Display color scale psscale -Coutcpt -D9/-1/18/0.2h -B2000:Topography: -K -O -V >> outps #Coastline and river database pscoast -R -JM -Dl -Ia/2/0/0/255 -A50 -W3 -K -O -V >>outps
#Psxy
psxy -JM -R -Sc0.1 -G255/0/0 -K -O -V <<END>> outps 116.28 39.54 121.26 31.12 121.31 25.02 END #TEXT pstext -R -JM -G0/0/0 -N -O -V <<END>>outps
104 63 35 0 10 CM China Topography
118 40.5 12 0 10 CM Beijing
126 31.5 12 0 10 CM Shanghai
123 26 12 0 10 CM Taipei
END
###########################################################
gv $outps &
20、Images clipped by coastlines
#!/bin/csh
######Author: didida 2006/11/22
######Purpose: Illustrates Clipping of images using coastline
set ingrd1 : China.grd
set ingrd2 : China_intens.grd
set outgrd1 : po1.grd
set outgrd2 : po2.grd
set outcpt1 : po1.cpt
set outcpt2 : po2.cpt
set outps : clipped2-1.ps
set range : 72/136/17.5/60
set proj : M6i
##########################################################
grdsample TOPO5.GRD -R72/136/17.5/60 -I10m -G$ingrd1 -V
makecpt -Csealand -T0/8000/100 -Z > $outcpt1
psbasemap -Rrange−JM6i−B5f5−X1−Y3−K−V>outps
grdgradient ingrd1−Gingrd2 -Nt1 -A120 -V
grdimage ingrd1−Iingrd2 -R -JM6i -Coutcpt1−K−O−V>>outps
psscale -Coutcpt1−D6.3/−1.5/13/0.3h−B1000−K−O−V>>outps
pscoast -R -JM -Dl -Sc -K -O -V >> $outps
echo "-10000 255 10000 150" > $outcpt2
grdgradient ingrd1−Goutgrd2 -Nt1 -A80 -V
grdimage ingrd1−Ioutgrd2 -JM -Coutcpt2−K−O−V>>outps
pscoast -R -JM -Q -B5f5:."Clipping of Images": -O -V >> $outps
pstext -R -JM -G255/0/0 -N -K -O -V << END >> $outps
122 25 15 0 15 CM Taipei
120 40 15 0 15 CM Beijing
125 31 15 0 15 CM Shanghai
END
########################################################
gv clipped2-1.ps &
21、Volume and spatial selections
#! /bin/csh
#####Author : didida 2006 11 30
#####Purpose : Volumes and Spatial Selection
#######################################################
set worldgrd : TOPO5.GRD
set ingrd1 : China.grd
set ingrd2 : China_intens.grd
set outgrd1 : mask.grd
set outgrd2 : mask_G.grd
set outgrd3 : tmp.grd
set output1 : centerpoint.d
set output2 : center.awk
set output3 : center.d
set outcpt : China18.cpt
set outps : China18.ps
set range : 72/136/17.5/60
set proj : M6i
set frame = a10f5/a5f5WSEn
########################################################
###Make Grdsample
makecpt -Chaxby -T-8000/8000/1000 -Z -V > outcptgrdsampleworldgrd -Gingrd1−Rrange -V -I10m
grdgradient ingrd1−Gingrd2−A0−M−Nt0.75b−Vgrdimageingrd1 -Iingrd2−Jproj -Rrange−Bframe -Coutcpt−X1i−Y1i−K−V>outps
psscale -Coutcpt−D8/−1/15/0.2h−B2000:Topography:−K−O−V>>outps
pscoast -R -Jproj−A2−W5−K−O−V>>outps
###Define Centerpoint Location of China
echo "90 33" > $output1
###draw a circle of radius = 1000 km centered on China
AWK′print$1,$2,12,0,1,"LB","Plateau"′output1 | pstext -R -JM -Gblue -N -O -K -D0.1i/0.1i >> outpsAWK '{print 1,2, 0, 1000, 1000}' output1|psxy−R−JM−O−K−SE−W0.5p/255/0/0−V>>outps
###Then draw 1000m contours and overlay 1000 m contour in green
grdcontour ingrd1−JM−R−C1000−Bframe -O -K -U/0.5i/-1i/"homework18" >> outpsgrdcontouringrd1 -JM -R -Bframe−C500−L990/1010−O−K−Dsm−Wc0.75p/0/255/0>>outps
pscoast -R -JM -O -K -W0.25p >> outpsAWK '{print 1,2, 0, 1000, 1000}' output1|psxy−R−JM−O−K−SE−W1p/255/0/0−V>>outps
#\rm -f sm_*.xyz # Only consider closed contours#
###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 output3foreachfile(sm∗.xyz)AWK -f output2file > $output3
end
###Only plot the ones within 5000km
gmtselect -R -JM -C1000/output1output3 >
psxy
-R -JM -O -K -SC0.04i -Gred -W0.5p >> outpspsxy−R−JM−O−K−SD0.1i−Gyellow−W0.5poutput1 >> $outps
###Then report the volume and area of these seamounts only by masking out data outside the
###1000 km-radius circle and then evaluate area/volume for the 1000m contour
grdmath -R -I10m 90 33 GDIST = outgrd1grdclipoutgrd1 -Sa1000/NaN -Sb1000/1 -Goutgrd2grdmathingrd1 outgrd2MUL=outgrd3
set area = `grdvolume outgrd3−C1000−Sk|cut−f3`setvolume=`grdvolumeoutgrd3 -C1000 -Sk | cut -f2`
pstext -R -JM -Gred -O -N << EOF >> outps11050901LMArea:area km@+2@+
110 45 9 0 1 LM Volumes: $volume km@+3@+
EOF
#########################################################
\rm -f $$ sm_*.xyz tmp.grd mask.grd center*
gv $outps &
22、Time-series along track
#! /bin/csh
###### HOMEWORK 3
###### Author: Dean 12/06/2006
###### Purpose: Plotting Time-Series along Tracks
#######################################################################################
#####set#####
set ingrd1 = TOPO5.grd
set ingrd2 = AtlanticOcean.grd
set ingrd3 = AtlanticOcean_intens.grd
set outps = AtlanticOcean.ps
set outcpt = Atlantic.cpt
set range = -80/20/0/60
set proj = M6i
set frame = a20f10/a20f10g10WSne
########################################################################################
gmtset ANOT_FONT_SIZE 10 LABEL_FONT_SIZE 10
makecpt -Crainbow -T-7000/3000/100 -Z -V > outcptgrdsampleingrd1 -Gingrd2−Rrange -V -I10m
grdgradient ingrd2−Gingrd3 -A0 -M -Nt0.75b -V
grdimage ingrd2−Iingrd3 -Rrange−Jproj -Bframe−Coutcpt -Y2i -P -K -V > outpspsscale−Coutcpt -D3i/-0.5i/6i/0.1ih -B1000:Topography: -K -O -V >> outpspscoast−R−Jproj -Dh -Ia/2/0/0/255 -A2 -W2 -K -O -V >> outps ######################################################################################## pswiggle track_d.*.xys -Rrange -Jproj−Bframe -Gpink -Z200 -W0.5p -S10/64/100/@~m@~rad \
-P -O >> outps ######################################################################################## gvoutps&
參考網站:http://topex.ucsd.edu/cgi-bin/get_data.cgi
獲得track_d*.xys檔--------------------------Magnetic Anomaly data
23、Gridding With Missing Data : Atlantic Ocean
#! /bin/csh
###### HOMEWORK 1
###### Author: Dean 11/17/2006
###### Purpose: Gridding and clipping when data are missing in Atlantic Ocean
#############################################################################
#####set#####
set inxyz = atlantic.xyz
set outgrd1 = AtlanticOcean.grd
set outgrd2 = ship_clipped10.grd
set outps = AtlanticOcean.ps
set output1 = AtlanticOcean.b
set output2 = AtlanticOcean10m.b
set range1 = -80/0/-89/89
set proj1 = M3
set frame = a40f10/a20f20WSne
########################################################################################
#####gmtconvert#####
gmtconvert inxyz−bo−V>output1
#####region#####
set region = `minmax output1 -I1 -bi3` #####nearneighbor##### nearneighboroutput1 -bi3 -Goutgrd1region -I10m -S200k -V
#####grdinfo#####
set info = `grdinfo -C -M outgrd1` #####grdcontour##### grdcontouroutgrd1 -Jproj1−Bframe -C1000 -A2000 -G2i -U -K -V > outps #####blockmedian##### blockmedianregion -bi3 output1−bo−I10m−V>output2
#####surface#####
surface output2−bi3−Goutgrd1 region -I10m -V #####psmask##### psmaskoutput1 -bi3 region−JM−I10m−T−G220−X2i−O−K−V>>outps
#####grdcontour#####
grdcontour outgrd1−JM−Bframe -C1000 -L-8000/0 -A2000 -G2i -O -K -V >> outps #####psmask##### psmaskoutput2 -bi3 region−JM−Bframe -I10m -X2i -O -K -V >> outps #####grdcontour##### grdcontouroutgrd1 -JM -C1000 -A2000 -L-8000/0 -G2i -O -K -V >> outps #####psmask##### psmask -C -O -K -V >>outps
#####grdclip#####
grdclip outgrd1−Sa−1/NaN−Goutgrd2 -V
#####grdcontour#####
grdcontour outgrd2−JM−Bframe -C1000 -A2000 -L-8000/0 -G2i -X2i -O -K -V >> outps #####pscoast##### pscoastregion -JM -G150 -W0.25p -O -K -V >> outps #####psxy##### echoinfo info| psxy−R−JM−Sa0.15i−W1p−O−K−V>>outps
#####pstext#####
echo "-2.5 6.8 25 0 1 CB Gridding with missing data:Atlantic Ocean" | \
pstext -R0/3/0/4 -Jx1i -N -O -V >> outps ######################################################################################## gvoutps &
24、Volume and spatial selections
#! /bin/csh
###### HOMEWORK 2
###### Author: Dean 11/27/2006
###### Purpose: Volumes and Spatial Selections in Atlantic Ocean
#######################################################################################
#####set#####
set ingrd1 = atlantic.grd
set ingrd2 = atlantic-intens.grd
set outgrd1 = maskatlantic.grd
set outgrd2 = maskatlantic-intens.grd
set outgrd3 = AtlanticOcean.grd
set outps = AtlanticOcean.ps
set output1 = AtlanticOcean.b
set output2 = AtlanticCenter.awk
set output3 = AtlanticCenter.b
set outcpt = Atlantic.cpt
set range = -80/10/-70/80
set proj = M3i
set frame1 = a20f10/a20f10g10WSne
set frame2 = a20f10/a20f10g10WSne
########################################################################################
#####Define location of South mid-ocean ridge
echo "-13 -30 " > output1 #####First generate gravity image w/ shading,label South mid-ocean ridge, #####and draw a ellipse of a=3000Km;b=1000Km centered on South mid-ocean ridge. makecpt -Crainbow -T-10000/6000/200 -Z >outcpt
grdgradient ingrd1−Nt1−A45−Gingrd2
grdimage ingrd1−Iingrd2 -Jproj−Coutcpt -Rrange−Bframe1 -P -K -Y2i -V > outpspsscale−Coutcpt -D1.5i/-0.5i/3i/0.1ih -B3000:Topography: -K -O -V >> outpspscoast−R−Jproj -A2 -W1 -K -O -V >> outpsAWK '{print 1,2, 12, 0, 1, "CM", "South mid-ocean ridge"}' output1|pstext−R−J −Gblue−N−O−K−D0.1i/0.1i>>outps
AWK′print$1,$2,0,3000,1000′output1 | psxy -R -Jproj−O−K−SE −W1p/255/0/0−V>>outps
##### Then draw -3000m contours and overlay -3000 m contour in green
grdcontour ingrd1−J−R−C3000−Bframe2 -O -K -U/0.5i/-1i -X4i >> outpsgrdcontouringrd1 -J -R -B -C500 -L-3100/-2900 -O -K -Dsm -Wc0.75p/0/255/0 >> outpspscoast−R−J−O−K−W0.25p>>outps
AWK′print$1,$2,0,3000,1000′output1 | psxy -R -Jproj−O−K−SE −W1p/255/0/0−V>>outps
#\rm -f sm_*.xyz # Only consider closed contours#
#########################################################################################
##### Now determine centers of each enclosed seamount > -3000m but only plot the ones
##### within South mid-ocean ridge
##### 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 -foutput3
foreach file (sm_*.xyz)
AWK−foutput2 file>>output3
end
#####Plot
gmtselect -R -Jproj−C3000/output1 output3>psxy−R−J−O−K−SC0.04i−Gred−W0.25p>>outps
psxy -R -J -O -K -SD0.1i -Gyellow -W0.25p output1>>outps
##### Then report the volume and area of these seamounts only by masking out data outside the
##### a ellipse of a=3000Km;b=1000Km centered on South mid-ocean ridge and then evaluate area/volume
#####for the -3000m contour
grdmath -R -I10m 47 24 GDIST = outgrd1grdclipoutgrd1 -Sa1500/NaN -Sb1500/1 -Goutgrd2grdmathingrd1 outgrd2MUL=outgrd3
set area = `grdvolume outgrd3−C1000−Sk|cut−f3`setvolume=`grdvolumeoutgrd3 -C1000 -Sk | cut -f2`
pstext -R -J -Gblue -O -N << EOF >> outps−55−601201LMAreas:area km@+2@+
-55 -62 12 0 1 LM Volumes: volume km@+3@+ EOF ################################################################################################## \rm -fsm∗.xyztmp.grdmask.grdcenter∗gvoutps&
页:
[1]