Matlab处理气象数据（十四）分段的逐点变化的DTR

% 两套数据逐点DTR，2000年前后对比
% 导入Tmax1,72*128*396
% 导入Tmax2,72*128*396
% 导入Tmin1,72*128*396
% 导入Tmin2,72*128*396
 
Ymax1=reshape(Tmax1,[72,128,12,33]);
Ymax1=nanmean(Ymax1,3);
Ymax1=squeeze(Ymax1); %<72x128x33 double>
 
Ymax2=reshape(Tmax2,[72,128,12,33]);
Ymax2=nanmean(Ymax2,3);
Ymax2=squeeze(Ymax2); 
 
Ymin1=reshape(Tmin1,[72,128,12,33]);
Ymin1=nanmean(Ymin1,3);
Ymin1=squeeze(Ymin1); 
 
Ymin2=reshape(Tmin2,[72,128,12,33]);
Ymin2=nanmean(Ymin2,3);
Ymin2=squeeze(Ymin2); 
 
DTR1=Ymax1-Ymin1;
DTR2=Ymax2-Ymin2;
 
DTR1=DTR1(:,:,3:33);
DTR1a=DTR1(:,:,1:19);
DTR1b=DTR1(:,:,20:31);
 
DTR2=DTR2(:,:,3:33);
DTR2a=DTR2(:,:,1:19);
DTR2b=DTR2(:,:,20:31);

%画NCEP数据2000年前的
for j=1:128
for i=1:72
y=DTR1a(i,j,:);
n = 19;
dt=1;
% 计算统计量
s = 0;
for p = 1:n-1
    for q = p+1:n
        s = s + sign( y(q) - y(p) );
    end
end
% 方差( assuming no tied groups )
v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;
% 检验统计量
if s == 0
z = 0;
elseif s > 0,
z = ( s - 1 ) / sqrt( v );
else
z = ( s + 1 ) / sqrt( v );
end
% should calculate Normal value here
nor = 1.64;
if abs( z ) < nor,
sya(i,j) = 0;
else
% 计算斜率
ndash = n * ( n - 1 ) / 2;
s = zeros( ndash,1 );
s=nan;
r = 1;
for p = 1:n-1
for q = p+1:n
s(r) = ( y(q) - y(p) ) / ( q - p ) / dt;
r = r + 1;
end
end
sDTR1a(i,j) = median( s );
end
end
end 
[x y]=meshgrid(72:0.5:135.5,18:0.5:53.5);
pcolor(x,y,sDTR1a);
shading flat;
load mycolor;
colormap(mycolor);
hold on
map=shaperead('E:\数据\边界\china_map1\maps\bou2_4l.shp');%加载省界带南海的边界线
bou2_4lx=[map(:).X];%提取经度
bou2_4ly=[map(:).Y];%提取纬度
provence=[bou2_4lx',bou2_4ly']; 
plot(bou2_4lx,bou2_4ly,'-k','LineWidth',1.2);%绘国界
axis([72 136 18 54]);%设置显示的经纬度范围
hold off
%最后调整区间为-0.2~0.2

%画NCEP数据2000年后的
for j=1:128
for i=1:72
y=DTR1b(i,j,:);
n = 12;
dt=1;
% 计算统计量
s = 0;
for p = 1:n-1
    for q = p+1:n
        s = s + sign( y(q) - y(p) );
    end
end
% 方差( assuming no tied groups )
v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;
% 检验统计量
if s == 0
z = 0;
elseif s > 0,
z = ( s - 1 ) / sqrt( v );
else
z = ( s + 1 ) / sqrt( v );
end
% should calculate Normal value here
nor = 1.64;
 
if abs( z ) < nor,
sya(i,j) = 0;
else
% 计算斜率
ndash = n * ( n - 1 ) / 2;
s = zeros( ndash,1 );
s=nan;
r = 1;
for p = 1:n-1
for q = p+1:n
s(r) = ( y(q) - y(p) ) / ( q - p ) / dt;
r = r + 1;
end
end
sDTR1b(i,j) = median( s );
end
end
end 
[x y]=meshgrid(72:0.5:135.5,18:0.5:53.5);
pcolor(x,y,sDTR1b);
shading flat;
load mycolor;
colormap(mycolor);
hold on
map=shaperead('E:\数据\边界\china_map1\maps\bou2_4l.shp');%加载省界带南海的边界线
bou2_4lx=[map(:).X];%提取经度
bou2_4ly=[map(:).Y];%提取纬度
provence=[bou2_4lx',bou2_4ly']; 
plot(bou2_4lx,bou2_4ly,'-k','LineWidth',1.2);%绘国界
axis([72 136 18 54]);%设置显示的经纬度范围
hold off
%最后调整区间为-0.2~0.2

%画观测数据2000年前的
for j=1:128
for i=1:72
y=DTR2a(i,j,:);
n = 19;
dt=1;
% 计算统计量
s = 0;
for p = 1:n-1
    for q = p+1:n
        s = s + sign( y(q) - y(p) );
    end
end
% 方差( assuming no tied groups )
v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;
% 检验统计量
if s == 0
z = 0;
elseif s > 0,
z = ( s - 1 ) / sqrt( v );
else
z = ( s + 1 ) / sqrt( v );
end
% should calculate Normal value here
nor = 1.64;
if abs( z ) < nor,
sya(i,j) = 0;
else
% 计算斜率
ndash = n * ( n - 1 ) / 2;
s = zeros( ndash,1 );
s=nan;
r = 1;
for p = 1:n-1
for q = p+1:n
s(r) = ( y(q) - y(p) ) / ( q - p ) / dt;
r = r + 1;
end
end
sDTR2a(i,j) = median( s );
end
end
end
[x y]=meshgrid(72:0.5:135.5,18:0.5:53.5);
pcolor(x,y,sDTR2a);
shading flat;
load mycolor;
colormap(mycolor);
hold on
map=shaperead('E:\数据\边界\china_map1\maps\bou2_4l.shp');%加载省界带南海的边界线
bou2_4lx=[map(:).X];%提取经度
bou2_4ly=[map(:).Y];%提取纬度
provence=[bou2_4lx',bou2_4ly']; 
plot(bou2_4lx,bou2_4ly,'-k','LineWidth',1.2);%绘国界
axis([72 136 18 54]);%设置显示的经纬度范围
hold off
%最后调整区间为-0.2~0.2

%画观测数据2000年后的
for j=1:128
for i=1:72
y=DTR2b(i,j,:);
n = 12;
dt=1;
% 计算统计量
s = 0;
for p = 1:n-1
    for q = p+1:n
        s = s + sign( y(q) - y(p) );
    end
end
% 方差( assuming no tied groups )
v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;
% 检验统计量
if s == 0
z = 0;
elseif s > 0,
z = ( s - 1 ) / sqrt( v );
else
z = ( s + 1 ) / sqrt( v );
end
% should calculate Normal value here
nor = 1.64;
if abs( z ) < nor,
sya(i,j) = 0;
else
% 计算斜率
ndash = n * ( n - 1 ) / 2;
s = zeros( ndash,1 );
s=nan;
r = 1;
for p = 1:n-1
for q = p+1:n
s(r) = ( y(q) - y(p) ) / ( q - p ) / dt;
r = r + 1;
end
end
sDTR2b(i,j) = median( s );
end
end
end
[x y]=meshgrid(72:0.5:135.5,18:0.5:53.5);
pcolor(x,y,sDTR2b);
shading flat;
load mycolor;
colormap(mycolor);
hold on
map=shaperead('E:\数据\边界\china_map1\maps\bou2_4l.shp');%加载省界带南海的边界线
bou2_4lx=[map(:).X];%提取经度
bou2_4ly=[map(:).Y];%提取纬度
provence=[bou2_4lx',bou2_4ly']; 
plot(bou2_4lx,bou2_4ly,'-k','LineWidth',1.2);%绘国界
axis([72 136 18 54]);%设置显示的经纬度范围
hold off
%最后调整区间为-0.2~0.2

得到如下结果：

NCEP数据1981-1999年（左）和2000-2011年（右）的DTR逐点趋势分析

观测数据1981-1999年（左）和2000-2011年（右）的DTR逐点趋势分析

---

相关链接：
Matlab处理气象数据——目录