毕业季要到了，毕业论文也要到了！希望能帮助到你

这个程序是做相关场的。这个的相关场是年数据的，也就是经过计算有N年的夏季数据和N年的nino某一季节的数据！

注意这里的NC文件层数是只有一层的，如果涉及多层要在NC文件读取那里改，一般NC文件都是从[ 经度  纬度 层数 时间 ] 这种顺序摆放的。

part1 输入文件的相关信息

这部分主要是输出相关NC文件和nino指数的信息，具体说明在注释里面

注意这里面的数据需要时间连续的数据

在这里就是基本上要改的内容了，都放在这里，还有不同NC文件的格式也要改，所以在NC文件读取那里就要改改

dt是滞后时间，改变dt可以做滞后相关(NC文件滞后dt年)

%% 相关场%相关说明%相关场%适用于全球数据，如果是海洋数据可以用patch覆盖陆地%第80行的m_coast可以加入('patch',[0.8 0.8 0.8])这个属性填充陆地%适用于只有一层的数据，如果多层数据要在NC文件读取那里改%t检验过检标记在最后面，如果没有需要可以直接删除%作者：施根号%% 输入的相关数据clc;clear;close all;dt = 0;                               %滞后年数reliability = 0.1;                    %信度检验，这里可以0.1，0.05或者0.01boundary = [0 360 -90 90];            %范围[起始经度 终止经度 起始纬度 终止纬度]time_start = datenum('01-01-1979');   % 月-日-年，如果是月数据再日期都是01time_end   = datenum('12-01-2018');   % 月-日-年，记得不要超出NC文件的时间inpath  = 'H:\NC文件.nc';  %NC文件位置inpath2 = 'H:\excel文件.xlsx';     %SSN指数文件位置，excel格式，记得只包含数据inpath3 = 'H:\图片输出位置';        %图片储存位置varname = 'NC文件的变量名';         %NC文件的变量名title_plot = 'XX-XX相关';          %标题及文件名

part2 数据的读取

2.1 nino指数的读取

因为大家的文件读写习惯不同，所以这里的excel文件是只有数据，没有其他说明的，比如这样：

也可以用其他格式的文件，对应换换函数就OK

%% nino指数的读取nino = xlsread(inpath2);len = length(nino);  %读取数据的年数

2.2 NC文件数据的读取

time_nc = time_nc/24 + datenum('1-jan-1800');%或者：time_nc = time_nc + datenum('1-jan-1800');

这里要特别注意这一句，matlab的时间戳和NC文件的时间戳是不一样的，所以要转换。我们可以ncdisp('文件路径')查看时间信息，如下图：

图1 units = 'hours since 1800-01-01 00:00:0.0'

图2 units= 'days since 1800-1-1 00:00:00'

第一，一些NC文件的时间系统是小时制的，也就是从说1800-1月-1日开始算到现在有多少小时，所以要除以24转换为天数；

就像图1，如果这个units = 'hours since 1800-01-01 00:00:0.0'这里是hours那就是要除以24

就像图2，如果是units= 'days since 1800-1-1 00:00:00'，这里面是days，就是以天数计算，所以不用除以24就可以

第二是matlab是从公元0000年-01月-01日开始计算到现在多少天的，而NC文件是从1800年1月1日算起，所以要加上datenum('1-jan-1800')以转换为matlab的时间戳。注意一些NC文件是从从1900年1月1日算起，我们可以从units = 'hours since 1800-01-01 00:00:0.0'这个获取相关信息

%% ---读取NC文件数据--ncid=netcdf.open(inpath,'NC_NOWRITE'); %打开nc文件lon = ncread(inpath,'lon');            %查阅经度信息loncount = length(lon);                %查阅经度的精度(有多少格点)lat = ncread(inpath,'lat');            %查阅纬度信息latcount = length(lat);                %查阅纬度精度(有多少格点)%查找绘制范围对应的所在矩阵的位置，相当于截取一小段矩阵lon_scope =  find(lon >= boundary(1) & lon<=boundary(2));lat_scope =  find(lat >= boundary(3) & lat<=boundary(4));%绘制范围的数据量，重新定义边界lon_number = length(lon_scope);lat_number = length(lat_scope);lat_1=lat(lat_scope(end):-1:lat_scope(1));lon_1=lon(lon_scope(1):lon_scope(end)); [plon,plat]=meshgrid(lon_1,lat_1);%重新定义边界boundary = [double(lon_1(1)) double(lon_1(end)) double(lat_1(1)) double(lat_1(end))]; %--------读取--------time_nc = ncread(inpath,'time');time_nc = time_nc/24 + datenum('1-jan-1800');time_start = find(time_nc == time_start);time_end = find(time_nc==time_end);time_number = length(time_nc(time_start : time_end));start = [lon_scope(1),lat_scope(1),time_start];    %初始位置count = [lon_number,lat_number,time_number];        %读取范围 stride1 = [1,1,1];                        %读取步长sst(:,:,:) = ncread(inpath,varname ,start,count,stride1);netcdf.close(ncid);      %关闭nc文件

part3 相关系数的计算

1. 里面的矩阵t是检验其过不过检验的指标。例如：t<0.1则通过了显著性水平α=0.1，在代码的开头有reliability = 0.1就是这个信度指标，可以选择0.01，0.05或者0.01。

2. 经过计算后的NC文件必须和nino指数的时间层数要对应上，比如有40年nino指数，那么经过计算后NC文件的时间层数也必须要有40年，若是要其他季节的平均可以在这次基础上改循环

%% 计算%----夏季平均-----for i = 1:len    sst_w(:,:,i) = (sst(:,:,6+ 12*(i-1)) + sst(:,:,7+ 12*(i-1)) + sst(:,:,8+ 12*(i-1)))/3;endsst_1 = sst_w(:,:,1+dt:end);nino_1 = nino(1:end-dt);n = length(nino_1);% ----计算相关系数及T检验---for i = 1:lon_number    for j = 1:lat_number        s = reshape(sst_1(i,j,:),n,1);        [r_temp,t_temp] = corrcoef(nino_1,s);        r(i,j) = r_temp(1,2);  %相关系数        t(i,j) = t_temp(1,2);  %检验    endendr =  rot90(r);t =  rot90(t);

part4 绘图

这里我们要用到m_map，如果还没装m_map可以先安装

如果是海洋数据(比如海表面温度SST)，那么第6行的m_coast可以加入('patch',[0.8 0.8 0.8])这个属性填充陆地

然后下面代码的第15行被注释的是把过了90%检验的地方圈起来，如果有需要可以取消注释

%% 绘图figure()m_proj('Equidistant','lat',[boundary(3) boundary(4)],'lon',[boundary(1) boundary(2)]);m_pcolor(plon,plat,r);              %添加我们要画的内容shading interp;m_coast('linewidth',1,'color','k');  %绘制海岸线% 添加边框m_grid('xtick',6,'ytick',5,'tickdir','in','yaxislocation','left','fontsize',10);title(title_plot,'fontsize',15); %标题colormap(jet);caxis([-1,1]);colorbar('location','Eastoutside');hold on%----------t检验标记，若不需要则删除-------------------------------------%[cs,h]=m_contour(plon,plat,t,[10,reliability],'r','linewidth',1.5);for i = 1:2:lon_number    for j = 1:2:lat_number        if t(j,i) <= reliability            m_plot(plon(1,i),plat(j,1),'k','marker','.','MarkerSize',2)            hold on;        end    endend%--------------删除到这------------------------------------------------cd (inpath3)                   %打开图片储存位置saveas(gcf,title_plot, 'png')  %保存图片

输出结果：

笔者水平有限，若有疏漏恳请指出！

如果对您有帮助可以点下文末的赞赏吗谢谢各位读者~

-END-

附：完整代码

%% 相关场%相关说明%相关场%适用于全球数据，如果是海洋数据可以用patch覆盖陆地%第80行的m_coast可以加入('patch',[0.8 0.8 0.8])这个属性填充陆地%适用于只有一层的数据，如果多层数据要在NC文件读取那里改%t检验过检标记在最后面，如果没有需要可以直接删除%作者：施根号%% 输入的相关数据clc;clear;close all;dt = 0;                               %滞后年数reliability = 0.1;                    %信度检验boundary = [0 360 -90 90];            %范围[起始经度 终止经度 起始纬度 终止纬度]time_start = datenum('01-01-1979');   % 月-日-年，如果是月数据再日期都是01time_end   = datenum('12-01-2018');   % 月-日-年，记得不要超出NC文件的时间inpath  = 'H:\NC文件.nc';  %NC文件位置inpath2 = 'H:\excel文件.xlsx';     %SSN指数文件位置，excel格式，记得只包含数据inpath3 = 'H:\图片输出位置';              %图片储存位置varname = 'NC文件的变量名';                      %NC文件的变量名title_plot = 'XX-XX相关';               %标题及文件名%% nino指数的读取nino = xlsread(inpath2);len = length(nino);  %读取数据的年数%% ---读取NC文件数据--ncid=netcdf.open(inpath,'NC_NOWRITE'); %打开nc文件lon = ncread(inpath,'lon');            %查阅经度信息loncount = length(lon);                %查阅经度的精度(有多少格点)lat = ncread(inpath,'lat');            %查阅纬度信息latcount = length(lat);                %查阅纬度精度(有多少格点)%查找绘制范围对应的所在矩阵的位置，相当于截取一小段矩阵lon_scope =  find(lon >= boundary(1) & lon<=boundary(2));lat_scope =  find(lat >= boundary(3) & lat<=boundary(4));%绘制范围的数据量，重新定义边界lon_number = length(lon_scope);lat_number = length(lat_scope);lat_1=lat(lat_scope(end):-1:lat_scope(1));lon_1=lon(lon_scope(1):lon_scope(end)); [plon,plat]=meshgrid(lon_1,lat_1);%重新定义边界boundary = [double(lon_1(1)) double(lon_1(end)) double(lat_1(1)) double(lat_1(end))];%--------读取--------time_nc = ncread(inpath,'time');time_nc = time_nc/24 + datenum('1-jan-1800');time_start = find(time_nc == time_start);time_end = find(time_nc==time_end);time_number = length(time_nc(time_start : time_end));start = [lon_scope(1),lat_scope(1),time_start];    %初始位置count = [lon_number,lat_number,time_number];        %读取范围 stride1 = [1,1,1];                        %读取步长sst(:,:,:) = ncread(inpath,varname ,start,count,stride1);netcdf.close(ncid);      %关闭nc文件%% 计算%----夏季平均-----for i = 1:len    sst_w(:,:,i) = (sst(:,:,6+ 12*(i-1)) + sst(:,:,7+ 12*(i-1)) + sst(:,:,8+ 12*(i-1)))/3;endsst_1 = sst_w(:,:,1+dt:end);nino_1 = nino(1:end-dt);n = length(nino_1);% ----计算相关系数及T检验---for i = 1:lon_number    for j = 1:lat_number        s = reshape(sst_1(i,j,:),n,1);        [r_temp,t_temp] = corrcoef(nino_1,s);        r(i,j) = r_temp(1,2);  %相关系数        t(i,j) = t_temp(1,2);  %检验    endendr =  rot90(r);t =  rot90(t);%% 绘图figure()m_proj('Equidistant','lat',[boundary(3) boundary(4)],'lon',[boundary(1) boundary(2)]);m_pcolor(plon,plat,r);              %添加我们要画的内容shading interp;m_coast('linewidth',1,'color','k');  %绘制海岸线% 添加边框m_grid('xtick',6,'ytick',5,'tickdir','in','yaxislocation','left','fontsize',10);title(title_plot,'fontsize',15); %标题colormap(jet);caxis([-1,1]);colorbar('location','Eastoutside');hold on%----------t检验标记，若不需要则删除-------------------------------------%[cs,h]=m_contour(plon,plat,t,[10,reliability],'r','linewidth',1.5);for i = 1:2:lon_number    for j = 1:2:lat_number        if t(j,i) <= reliability            m_plot(plon(1,i),plat(j,1),'k','marker','.','MarkerSize',2)            hold on;        end    endend%--------------删除到这------------------------------------------------cd (inpath3)                   %打开图片储存位置saveas(gcf,title_plot, 'png')  %保存图片

‍