matlab滴m_map学习心得4

%% 读数据
file='C:\Users\Administrator\Desktop\专题四 海洋水质监测\COMS_GOCI_L2A_GA_20200322031640_L2.he5';
ncdisp(file);
longitude=ncread(file,'/navigation_data/longitude');
latitude=ncread(file,'/navigation_data/latitude');
Rrs_412=ncread(file,'/geophysical_data/Rrs_412');
Rrs_443=ncread(file,'/geophysical_data/Rrs_443');
Rrs_490=ncread(file,'/geophysical_data/Rrs_490');
Rrs_555=ncread(file,'/geophysical_data/Rrs_555');
Rrs_660=ncread(file,'/geophysical_data/Rrs_660');
Rrs_680=ncread(file,'/geophysical_data/Rrs_680');
solz=ncread(file,'/geophysical_data/solz');
[LON,LAT] = meshgrid(lon_real,lat_real);
%% 设置三维变量空间
Rrs_670=(Rrs_660+Rrs_680)/2;
Rrs(:,:,1)=Rrs_412;
Rrs(:,:,2)=Rrs_443;
Rrs(:,:,3)=Rrs_490;
Rrs(:,:,4)=Rrs_555;
Rrs(:,:,5)=Rrs_660;
Rrs(:,:,6)=Rrs_670;
Rrs(:,:,7)=Rrs_680;
wavelength = [412,443,490,555,660,670,680];
aw = [0.004562,0.00707,0.0150,0.0596,0.41,0.439,0.465];
bw = [0.00666,0.00478,0.0031,0.00185,0.0008,0.0008,0.0007];
bbw = [0.00335,0.00245,0.00158,0.0009,0.00044,0.00041,0.00038];
%% 公式
for i=1:80
    for j=1:90
        Rrs1 = reshape(Rrs(i,j,:),[1,7]);
        rrs=Rrs1./(0.52+1.7*Rrs1);
        u=(-0.089+sqrt(0.089^2+4*0.125*rrs))/(2*0.125);
        ka=log10((rrs(2)+rrs(3))/(rrs(4)+(5*rrs(6)*rrs(6))/rrs(3)));
        yita=2.0*(1-1.2*exp(-0.9*rrs(2)/rrs(4)));
        if Rrs1(6)<0.0015
            a1=0.0596+10^(-1.146-1.366*ka-0.469*ka^2);
            Bbp1=u(4)*a1./(1-u(4))-0.0009;
            Bbp=Bbp1*(555./wavelength).^yita;
        else
            a2=0.439+0.39*((rrs(6)/(rrs(2)+rrs(3)))^1.14);
            Bbp2=u(6)*a2./(1-u(6))-0.00041;
            Bbp=Bbp2*(670./wavelength).^yita;
        end
        Bb=bbw+Bbp;
        a=(1-u).*(Bbp+bbw)./u;
        Kd=(1+0.005.*solz(i,j)).*a+4.259*(1-(0.265*bbw./Bb)).*(1-0.52.*exp(-10.8*a)).*Bb;
        [M,I]=min(Kd);%M：最小值；I：最小值所在位置
        Zsd(i,j)=(1/(2.5*M))*log(abs(0.14-Rrs1(I))/0.013);
    end
end
%% 画图
figure(1)  
m_proj('miller','lon',[120.09  120.53],'lat',[35.9 36.28]); 
m_pcolor(longitude,latitude,Zsd);
m_gshhs_h('color','k');
m_coast('patch',[.7 .7 .7]);
m_grid('fontsize',12,'fontname','Times New Roman');
colormap(flipud(jet(256)));
ax=m_contfbar(.97,[.15 .85],[0 4],[0:0.001:4],'edgecolor','none','endpiece','no','ticklength',[.01 0.05]);
set(ax,'ytick',[0 1 2 3 4],'yticklabel',{'0','1','2','3','4'});
title(ax,'Zsd(m)','color','k','fontweight','bold','fontname','Times New Roman')  ;
m_ruler([.05 .3],.95,4,'color','b','linewid',2,'tickdir','out','ticklen',0.01,'fontsize',8);
m_northarrow(120.5,36.24,.03,'type',4);