matlab HYCOM数据画温度断面图

下载地址：
http://www.hycom.org/dataserver/gofs-3pt1/analysis

%%
clear,clc
%读取温度
clear,clc
temperature=[];
datadir='路径名'; %指定批量数据所在的文件夹
filelist=dir([datadir,'*.nc']); %指定批量数据的类型
k=length(filelist);
for s=1:k
  filename=[datadir,filelist(s).name];
  ncid=netcdf.open(filename,'NC_NOWRITE');
  t = ncread(filename,'water_temp'); 
  temperature = cat(3,temperature,t);
  netcdf.close(ncid);   % 关闭文件
end;

depth=[0 2 4 6 8 10 12 15 20 25 30 35 40 45 50 60 70 80 90 100 125 150 200]';
[Lat Lon]=meshgrid(lat,lon);
idx=find(lon>=121&lon<=124);idy=find(lat==28);idz=find(depth<=200);
TT=squeeze(temperature(idx,idy,idz,:));%去除纬度
TT=nanmean(TT,3);%在时间维度上求非nan的平均
TX=lon(idx);TY=lat(idy);TZ=depth(idz);%设定经纬度和深度的范围
[TZ,TX]=meshgrid(TZ,TX);
pcolor(TX,TZ,TT);%画温度断面图
set(gca,'YDir','reverse');shading interp;
caxis([10 21]);colormap('jet');cb=colorbar;

a1=find(isnan(TT(1,:))==1,1,'first');
pz(1)=TZ(1,a1);
a2=find(isnan(TT(2,:))==1,1,'first');
pz(2)=TZ(1,a2);
for i=3:length(TX(:,1))
    if sum(isnan(TT(i,:))==1)~=0
        %nan数为1,寻找逻辑数1的位置,返回第1个索引
        a=find(isnan(TT(i,:))==1,1,'first');
        pz(i)=TZ(1,a-1);
    else
        pz(i)=max(TZ(1,1:end-1));
    end
end
%画海底
p=patch([TX(:,1)' max(TX(:,1)) min(TX(:,1))],[pz max(TZ(1,:)) max(TZ(1,:))],'k');
set(p,'facecolor',0.5*[1 1 1])
xlabel('longitude(^o)');ylabel('depth(m)');ylabel(cb,'temperature');
title('温度断面图')