GMT计算二维网格均值、标准差和趋势

下面介绍GMT中blockmean和grdmath的统计计算功能，包含均值、标准差和趋势。

案例：高度计大气湿延迟的中国周边海域分布

以中国周边海域的高度计湿延迟为例，计算湿延迟的空间均值、标准差和趋势。

大气湿延迟：由于大气水汽对无线电传播速度的衰减造成，雷达高度计数据中对雷达测距的一个改正项。

Code
; @' W5 o+ D/ W* p  h/ i
$ b6 ]) Y8 y; ~. X/ o0 ~* b1 v' W+ z        
; E' o8 Y$ H! g4 |: n6 q+ v

[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

" A: `2 D2 X  ?* j. a3 |2 _( Z
结果

大气湿延迟和大气水汽和液态云有关，均值呈现出南高北低的自然特征。在北纬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

* k7 \) s2 V1 }9 P( Q