clc
clear all

% Ocean colour data from http://seawifs.gsfc.nasa.gov/SEAWIFS.html
%
% Take a 4km weakly average dataset and plot a map for the Strait of
% Georgia and outer coast. Note that most of this code is used
% for reading in and subsetting the data.
path='D:\Contour test\';    %貼上資料及程式所放的資料夾路徑
cd(path);                           %前往路徑的資料夾

%製作資料夾
mkdir 140E-180E-15S-15N

%指定經緯度及路徑
l.path=[path '140E-180E-15S-15N\'];                                              %(1) 140E-180E-15S-15N 最大範圍
l.lon=[130 210];         l.lat=[-20 20];

b = dir('*L3m**NSST*');
for i=1:length(b);

    % Note - This is probably not the most efficient way to read and
    %        handle HDF data, but I don't usually do this...

        % First, get the attribute data
        PI=hdfinfo(b(i).name);
        % And write it into a structure
        pin=[];
        for k=1:60,
          nm=PI.Attributes(k).Name;nm(nm==' ')='_';
          if isstr(PI.Attributes(k).Value),
            pin=setfield(pin,nm,PI.Attributes(k).Value);
          else
            pin=setfield(pin,nm,double(PI.Attributes(k).Value));
          end
        end;  

        
    for path_j=1
                
        % lon/lat of grid corners
        lon=[pin.Westernmost_Longitude:pin.Longitude_Step:pin.Easternmost_Longitude];
        lat=[pin.Northernmost_Latitude:-pin.Latitude_Step:pin.Southernmost_Latitude];
        
        if l.lon(2)>=180
            lon(lon<0) = lon(lon<0)+360;
            lon = [lon(4321:8640) lon(1:4320)];
        end

        % Get the indices needed for the area of interest
        [mn,ilt]=min(abs(lat-max(l(path_j).lat)));
        [mn,ilg]=min(abs(lon-min(l(path_j).lon)));
        ltlm=ceil(diff(l(path_j).lat)/pin.Latitude_Step);
        lglm=ceil(diff(l(path_j).lon)/pin.Longitude_Step);

        % load the subset of data needed for the map limits given
        if l.lon(2)>=180
            P=hdfread(b(i).name,'l3m_data','Index',{[],[],[]});
            P=[P(:,4321:8640) P(:,1:4320)];
            P=P(ilt:(ilt+ltlm),ilg:(ilg+lglm));
        else
            P=hdfread(b(i).name,'l3m_data','Index',{[ilt ilg],[],[ltlm+1 lglm+1]});           %+1為避免周圍出現空白
        end
        
        % Convert data into log(Chla) using the equations given. Blank no-data.
        P=double(P);

        P(P==65535)=NaN;
        P=(pin.Slope*P+pin.Intercept);

        LT=lat(ilt+[0:ltlm]);LG=lon(ilg+[0:lglm]);
        [Plg,Plt]=meshgrid(LG,LT);

        % Draw the map...
        cd(l(path_j).path)
        
        clf;                                                                              %開啟畫圖程式
        m_proj('miller','lon',l(path_j).lon,'lat',l(path_j).lat);                                     %設定投影方式
        m_pcolor(Plg,Plt,P);                                                              %畫出圖形
        shading flat;                                                                     %去除網格
        hold on
        m_contourf(Plg,Plt,P,[29 29],'color','k');  %畫某一特定值等值線時，[ ]內兩數字要相同
        hold off
        m_gshhs_i('color','k');                                                           
        m_grid('tickdir','out');                                                          

% m_grid('linewi',2,'tickdir','out');
% m_grid('box','fancy','tickdir','out');
        set(gca,'Clim',[20,30]);                                            %限定上限30下限20

        h=colorbar;
        set(get(h,'ylabel'),'String','Temperature (^{o}C)');
        
        set(h,'ytick',[ 20:2:30 ],...
              'yticklabel',[ 20:2:30 ],...
              'tickdir','out','fontsize',9);
        
        title(['MODIS SST   ' datestr(datenum(pin.Period_Start_Year,1,0)+pin.Period_Start_Day,26) ' -> ' ...
                              datestr(datenum(pin.Period_End_Year,1,0)+pin.Period_End_Day,26)],...
               'fontsize',14,'fontweight','bold');
           
        saveas(gcf,...
           ['SST4-' datestr(datenum(pin.Period_Start_Year,1,0)+pin.Period_Start_Day,29) '--' ...
            datestr(datenum(pin.Period_End_Year,1,0)+pin.Period_End_Day,29) '.png'],'png');

    cd(path);
    end

    disp('done.');

end

        
% clear all