|
下面介绍GMT中blockmean和grdmath的统计计算功能，包含均值、标准差和趋势。 案例：高度计大气湿延迟的中国周边海域分布以中国周边海域的高度计湿延迟为例，计算湿延迟的空间均值、标准差和趋势。 大气湿延迟：由于大气水汽对无线电传播速度的衰减造成，雷达高度计数据中对雷达测距的一个改正项。
Code
+ ? s1 l, n8 Q# k6 S0 J G; ^/ ]1 v
& R' \$ S# z3 ^9 H1 c, U[Python] 纯文本查看 复制代码
REM GMT EXAMPLE xxxx
REM 2021-03-14 [email]leiyang@fio.org.cn[/email]
REM Purpose: Using blockmean, grdmath to calculate the STD abd trend
gmt gmtset FORMAT_GEO_MAP = dddF MAP_FRAME_WIDTH=2p
gmt gmtset FONT_ANNOT_PRIMARY 7p,Helvetica,black FONT_LABEL 7p,Helvetica,black
set ps=jason_wet_global_new.ps
set R=105/130/15/45
set J=M5c
REM Calculate mean value and plot mean wet path delay
gmt makecpt -Cvik -T50/350/5> t.cpt
REM use Jason-2 data of cycle 1-294 (6years), downloaded from RADS.
rem gawk "!/NaN/ && NR>13 && $3>100 {print $3,$2,(-$4)*1000}" D:\rads\wet\radsdata_ok\j2*c*.asc > model.dd
gmt blockmean model.dd -R%R% -I20m >tmp_.dd
gmt surface tmp_.dd -R%R% -I20m -Ggridded2_j2_s.nc -T0.25
REM gridded2_j2_s.nc is the averaged file and will be used in next step.
gmt grdimage gridded2_j2_s.nc -R%R% -J%J% -K -Ct.cpt -Y13c> %ps%
gmt pscoast -R -J -Dl -A10000/0/1 -Bag --FONT_TITLE=10p -Gwhite -W0.1p --MAP_ANNOT_OBLIQUE=45 -O -K >> %ps%
gmt psscale -DjBC+o0c/-1.4c+w1.5i/0.08i -R -J -Ct.cpt -Bxaf -By+lmm -I -O --FONT_ANNOT_PRIMARY=7p -K>> %ps%
echo 107 42 a)mean of WPD| gmt pstext -F+f8p,1,black+jLB -R -J -O -N -Gwhite -K >> %ps%
rem Calculate STD and PLOT
gmt makecpt -Cvik -T50/150/1> t.cpt
gmt blockmean model.dd -As -E -R -I20m -Ggridded2_j2_std.nc
gmt grd2xyz gridded2_j2_std.nc | gawk "!/NaN/ {print $1,$2,$3}" | gmt surface -R -I20m -Ggridded2_j2_std2.nc -T0.25
REM gridded2_j2_std2.nc is the STD file and will be used in next step.
gmt grdimage gridded2_j2_std2.nc -R -J -K -Ct.cpt -O -X7c >> %ps%
gmt pscoast -R -J -Dl -A10000/0/1 -Bag --FONT_TITLE=10p -Gwhite -W0.1p --MAP_ANNOT_OBLIQUE=45 -O -K >> %ps%
gmt psscale -DjBC+o0c/-1.4c+w1.5i/0.08i -R -J -Ct.cpt -Bxaf -By+lmm -I -O --FONT_ANNOT_PRIMARY=7p -K>> %ps%
echo 107 42 b)STD of WPD| gmt pstext -F+f8p,1,black+jLB -R -J -O -N -K -Gwhite >> %ps%
rem Calculate trend and PLOT
gmt makecpt -Cvik -T-30/30/1> t.cpt
gmt grdsample gridded2_j2_s.nc -I60m -Ggridded2_j2_s_re.nc
gmt grdmath gridded2_j2_s_re.nc DDY = slope2.nc
gmt grdmath gridded2_j2_s_re.nc DDX = slope3.nc
gmt grdmath gridded2_j2_s_re.nc D2DY2 = slope1.nc
gmt grdimage slope2.nc -R -J -K -Ct.cpt -O -X-7c -Y-10c>> %ps%
gmt pscoast -R -J -Dl -A10000/0/1 -Bag --FONT_TITLE=10p -Gwhite -W0.1p --MAP_ANNOT_OBLIQUE=45 -O -K >> %ps%
gmt psscale -DjBC+o0c/-1.4c+w1.5i/0.08i -R -J -Ct.cpt -Bxaf -By+lmm/@. -I -O --FONT_ANNOT_PRIMARY=7p -K>> %ps%
echo 107 42 c)trend of WPD| gmt pstext -F+f8p,1,black+jLB -R -J -O -N -Gwhite -K >> %ps%
gmt grdimage slope3.nc -R -J -K -Ct.cpt -O -X7c >> %ps%
gmt pscoast -R -J -Dl -A10000/0/1 -Bag --FONT_TITLE=10p -Gwhite -W0.1p --MAP_ANNOT_OBLIQUE=45 -O -K >> %ps%
gmt psscale -DjBC+o0c/-1.4c+w1.5i/0.08i -R -J -Ct.cpt -Bxaf -By+lmm/@. -I -O --FONT_ANNOT_PRIMARY=7p -K>> %ps%
echo 107 42 d)trend of WPD| gmt pstext -F+f8p,1,black+jLB -R -J -O -N -Gwhite >> %ps%
gmt psconvert %ps% -A -P -Tf
5 Q, X" d2 c" Y7 F1 f
/ p, |. ` P9 f4 M8 f8 [8 U$ K4 b6 x1 M7 y% Q- | E
结果
大气湿延迟和大气水汽和液态云有关，均值呈现出南高北低的自然特征。在北纬25-35度之间，时间变化（STD）最为剧烈，说明这段区间可能四季分明，不过度干燥，不过分潮湿，适合人类居住（xia che）。变化趋势图表示南北和东西方向一度距离上湿延迟的变化，因为网格的空间单位是度，所以这里一阶导数结果表示一度空间距离上的湿延迟变化量。 要点blockmean model.dd -As -E中的-As设置输出为STD。如果不做设置，默认输出的是均值。 gmt grdmath gridded2_j2_s_re.nc DDY = slope2.nc中的DDY或者DDX设置一个固定方向的一阶导数计算，因为网格被提前预设为1度，因此计算结果是1度距离上的变化率。grdmath也可以计算二阶导数，使用D2DY2或者D2DX2。 数据可以直接使用小编从RADS数据库中提取的文件model.dd，它是 Jason-2 cycle 1-294的全部模型湿延迟原始数据汇总： [Python] 纯文本查看 复制代码 model.dd：[url=https://www.jianguoyun.com/p/DWQl1jsQ2PCQBxjWsYwE]https://www.jianguoyun.com/p/DWQl1jsQ2PCQBxjWsYwE[/url]
当然也可以下载原始数据，然后使用下面一行代码提取数据： [Python] 纯文本查看 复制代码 nt $3,$2,(-$4)*1000}" D:\rads\j2*c*.asc > model.dd
3 ^6 {: h4 F) X: l9 t1 K( o" G | 
