0%

使用GMT计算二维网格数据的均值、标准差和趋势

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

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

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

code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
REM             GMT EXAMPLE xxxx
REM 2021-03-14 leiyang@fio.org.cn
REM Purpose: Using blockmean, grdmath to calculate the STD abd trend
REM
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

结果

大气湿延迟和大气水汽和液态云有关,均值呈现出南高北低的自然特征。在北纬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的全部模型湿延迟原始数据汇总:
model.dd: https://www.jianguoyun.com/p/DWQl1jsQ2PCQBxjWsYwE