Matlab处理气象数据(十三)逐点变化的DTR及其异常
程序员文章站
2022-07-14 11:54:11
...
% 两套数据逐点DTR
% 导入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;
画两套数据的DTR变化趋势图
%NCEP的
for j=1:128
for i=1:72
if isnan(DTR1)==1
continue
else
y=DTR1(i,j,:);
n = 33;
dt=1;
% 计算斜率
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
sDTR1(i,j) = median( s );
end
end
end
pcolor(sDTR1);
shading flat;
%添加边界线
[x y]=meshgrid(72:0.5:135.5,18:0.5:53.5);
pcolor(x,y,sDTR1);
shading flat;
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
%最后调整colormap颜色区间为-0.1~0.1
%观测数据的
for j=1:128
for i=1:72
if isnan(DTR2)==1
continue
else
y=DTR2(i,j,:);
n = 33;
dt=1;
% 计算斜率
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
sDTR2(i,j) = median( s );
end
end
end
pcolor(sDTR2);
shading flat;
%添加边界线
[x y]=meshgrid(72:0.5:135.5,18:0.5:53.5);
pcolor(x,y,sDTR2);
shading flat;
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.1~0.1
相关链接:
Matlab处理气象数据——目录