 % proc_peru_phos_sil2_3  Process the Peru data sent from Ruth, see Francisco, create a two panel figure: contours and sample locations.
%  The data files were initially processed by Dan Cosgrove.  Data are on \\tornado\diatom\peru\Ruth
%  Data were extracted from the ACCESS database peru_historic_interpolated.mdb
%  Use the chkdist variable to check the distribution of data.  You may want to modify the color axis and levels based upon this 
%  information.  
%  This program was from John Ryan and modified to produce the Peru plots.  09/2003
%
% rpm

    close all; clear all; clear global; pack;
    
    %load sil_1000_1;
    %load sil_1000_2; 
    %sil_m1 = medfilt1(sil_1000_1, 6);                           %po4_m1 = medfilt1(pho4_100_1, 5);    
    %sil_m2 = medfilt1(sil_1000_2, 6);  
    
    InputData=load('peru_contours2.asc'); %  load p; Contains Peru data for depths < 201
    InputDataPoints = size(InputData)
    SeasonData = InputData(:,1);
    
    SeasonId = SeasonData(:,1)-floor(SeasonData(:,1));
    SeasonIdPoints = size(SeasonId)
    
    lat = InputData(:,2); 
    lon = InputData(:,3); 
    
    chkdist=0;  % flag used to process the parameters to check the distribution, see line 20
    load jetplus3;
    load exclude;
   
    
    PlotNames={'Phosphate (uM)';'Silicate (uM)'};  % plot labels  
    depthchk=0;
    depthn= '1000';
   
    for seasonc =4:4,          %  seasons 1-s.spring, 2-s.summer, ...  5-annual
       
    %-------------------------------seasons -----------------
    season_chk=seasonc;       
    sname = {'Spring';'Summer';'Fall';'Winter'};
   
    % 0.73=< s.spring  <0.96                  exclude  (season(:)>=0.96 | season(:)<0.73) 
    % O.O =< s.summer  <0.22 & 0.96< s.summer <0.9999  (season(:)>=0.22 & season(:)<0.96)
    % 0.22=< s.fall    <0.47                           (season(:)>=0.47 | season(:)<0.22)
    % 0.47=< s. winter <0.73                           (season(:)>=0.73 | season(:)<0.47)                   
    if season_chk==1,
        idx=find(InputData(:,3)>-50 | InputData(:,2)>0 | InputData(:,4)<depthchk | InputData(:,4)>depthchk | SeasonId(:)>=0.96 | SeasonId(:)<0.73); 
        InputData(idx,:)=[];  % filter data; modify p(:,4) to change depth
        SeasonName = 'Spring'
    end 
    if season_chk==2,
        idx=find(InputData(:,3)>-50 | InputData(:,2)>0 |InputData(:,4)<depthchk | InputData(:,4)>depthchk |(SeasonId(:)>=0.22 & SeasonId(:)<0.96)); 
        InputData(idx,:)=[];  % filter data; modify p(:,4) to change depth
        SeasonName = 'Summer'
    end 
    if season_chk==3,
        idx=find(InputData(:,3)>-50 | InputData(:,2)>0 | InputData(:,4)<depthchk | InputData(:,4)>depthchk | SeasonId(:)>=0.47 | SeasonId(:)<0.22); 
        InputData(idx,:)=[];  % filter data; modify p(:,4) to change depth
        SeasonName = 'Fall'
    end    
    if season_chk==4,
        idx=find(InputData(:,3)>-50 | InputData(:,2)>0 | InputData(:,4)<depthchk | InputData(:,4)>depthchk | SeasonId(:)>=0.73 | SeasonId(:)<0.47); 
        InputData(idx,:)=[];  % filter data; modify p(:,4) to change depth
        SeasonName = 'Winter'
    end  
    if season_chk==5,
        idx=find(InputData(:,3)>-50 | InputData(:,2)>0 | InputData(:,4)<depthchk | InputData(:,4)>depthchk) ; 
        InputData(idx,:)=[];  % filter data; modify p(:,4) to change depth
        SeasonName = 'Annual'
    end  
    
    NoDataPoints = size(idx)
    time=InputData(:,1); 
    lat=InputData(:,2); 
    lon=InputData(:,3); 
    
    LattPos = lat;
    LongPos = lon;
    
    z=InputData(:,4); 
    pho=InputData(:,5); 
    sil=InputData(:,6);
    % filter outliers
    a=find(pho<=0); pho(a)=nan;
    size(a)
    a=find(sil<=0 ); sil(a)=nan;         %| no3>100); no3(a)=nan;
    size(a)
   
    paramss=str2mat('pho','sil');   % strings used in output file names, see line 119
    if chkdist 
        figure(1);clf; 
        subplot(1,3,1); hist(pho,50);
        subplot(1,3,2); hist(sil,50); return
    end
    %Color axis derived from distribution. Make modifications here for different depths.
    ColorScaleRange=[0.50 3.50; 3.00 96.0];        %(3 6 oxy, 0 10 nit) min and max for phos, silicate
    
    %CBT=log([0.4 .8 1.6 3.2]);                          % phos log
    %CBTL=['0.4';'0.8';'1.6';'3.2'];
    
    ColorBarTick2  = log([3.0 6.0 12.0 24.0 48.0 96.0]);                % sili log
    ColorBarLabel2 = ['3.00';'6.00';'12.0';'24.0';'48.0';'96.0'];  
    
     disp (' dor ')    
    
    % :::::::::::::::::: M A I N   L O O P :::::::::::::::::::::::
    % Cycle through the parameters using the variable 'count' to advance to the next column in the data file
    for count = 2:2,
        eval(['param=' paramss(count,:) ';']);
        nbin=100;
        rlon=linspace(min(lon),max(lon),nbin);
        rlat=[linspace(min(lat),max(lat),nbin)]';
        rparam=repmat(NaN,nbin-1); rparamN=rparam;
        for N=1:length(rlon)-1;
            for M=1:length(rlat)-1;
                idx=find(lon>=rlon(N) & lon<=rlon(N+1) & lat>=rlat(M) & lat<=rlat(M+1) & ~isnan(param));
                if ~isempty(idx); cparam=param(idx); rparam(M,N)=median(cparam); rparamN(M,N)=length(cparam); end
            end
        end     
        rparam=mf(rparam,3,1);
        rlon=rlon(1:length(rlon)-1); rlat=rlat(1:length(rlat)-1);
        rlon=rlon+median(diff(rlon)); rlat=rlat+median(diff(rlat));
        load exclude
        for N=1:9;
            eval(['X=x' int2str(N) '; Y=y' int2str(N) ';']); X=sort(X); Y=sort(Y);
            xidx=find(rlon>=X(1) & rlon<=X(2)); yidx=find(rlat>=Y(1) & rlat<=Y(2));
            rparam(yidx,xidx)=repmat(NaN,length(yidx),length(xidx));
        end
        [glon glat]=meshgrid(rlon,rlat);
        idx=find(finite(rparam));
        gparam=griddata(glon(idx),glat(idx),rparam(idx),rlon,rlat);
        gparam=mf(gparam,5,1); 
        tmp=min(gparam); minp = min(gparam(:));
        tmp=max(gparam); maxp = max(tmp);
        
        PlotAxesRange = [-84 -72 -18 -4]; 
        color_scale=ColorScaleRange(count,:);
       % if count==1                         %for phos only
        %    cax=log(cax); gparam=log(gparam);
        % end
        if count==2                         %for silcate only
           color_scale=log(color_scale); gparam=log(gparam);
       end
        
        
        FS=12; clear XTL YTL; deg=setstr(176); XT = [-84:2:-72]; YT = [-18:2:-4];
        for N=1:length(XT); XTL(N,:)=[int2str(abs(XT(N))) setstr(176) 'W']; end
        for N=1:length(YT); 
            if YT(N)>-10;
                YTL(N,:)=[' ' int2str(abs(YT(N))) setstr(176) 'S']; 
            else
                YTL(N,:)=[int2str(abs(YT(N))) setstr(176) 'S']; 
            end
        end
        %figure 1;
        Figure_h = eval(['figure(' int2str(count) ');']);           % Create Figure and get ist handle
        clf;                                                            % Clear Figure
        set(Figure_h,'units','pixels','position',[10 100 800 500],'color','w');
        colormap(jetplus3);
        % plot positions in Figure
        px=.07; py=.08; pw=.4; ph=.8;
        Position_1 = [px py pw ph];
        Position_2 = [px+pw+.07 py pw ph];
        BarPosition = [Position_1(1)+Position_1(3)-0.03, Position_1(2)+.2, .01, Position_1(4)-.3];
        %-----------------------------------------------------
        
        Plot1_Axes_h = axes('position', Position_1);
        m_proj('Equidistant Cylindrical', 'lon', PlotAxesRange([1 2]),'lat', PlotAxesRange([3 4])); 
        Contour_Lines = 10;
        [cc,Contour_h] = m_contourf(rlon,rlat,gparam, Contour_Lines);        
        caxis(color_scale);                                     % Scale color range with color_scale
        %set(Contour_Patch_h(:),'edgecolor',[0.7,0.7,0.7]);  
        
        %-----   Extract Contour Lines
        ss = size(Contour_h);
        Contours = ss(1)
        for ic = 1:Contours
            lx = get(Contour_h(ic), 'XData')
            cxline(1:size(lx), ic) = lx(:);
            ly = get(Contour_h(ic), 'YData');
            cyline(1:size(ly), ic) = ly(:); 
            ss = size(ly);
            line_points(ic) = ss(1);;
        end;
        %-----   Save the land patches once, and reuse later
        %m_gshhs_h('save', 'ZK_Coastline');                     % Save into file
        LandColor = [0.7, 0.7, 0.7];
        m_usercoast('ZK_Coastline', 'patch', LandColor)        % Read from file           
        SaveAxis = axis;            % save for the next plot
        
        %-----   here define (or read) the "NoDaTA" Patch
        xpatch = [-0.08, 0.08, 0.08, -0.08];
        ypatch = [-0.15, -0.25, -0.3, -0.25];
        PatchColor = [0.9, 0.9, 0.9];
        patch(xpatch, ypatch, PatchColor, 'edgecolor', PatchColor);
        
        %-----   Plot contour with selective colors
        for ic = 1:Contours
            llx = cxline(1:line_points(ic), ic);
            lly = cyline(1:line_points(ic), ic);
            IN = inpolygon(llx, lly, xpatch, ypatch);
            for i = 1:size(llx) - 1; 
                if IN(i) == 1 color = PatchColor;  else color = LandColor; end;
                line(llx(i:i+1), lly(i:i+1), 'Color', color);
            end;
        end;    
        set(Plot1_Axes_h,'Xtick',[],'Ytick',[]);                
        g=get(Plot1_Axes_h,'position');   
        Plot1_Axes_Over_h = axes('position',g); 
        axis(Plot1_Axes_Over_h, PlotAxesRange); 
        set(Plot1_Axes_Over_h,'color','none','visible','off');
        gry=[.7 .7 .7]; 
        yllocright=0; 
        lblaxes=1; 
        lcol='k'; 
        scfac=.7;
        mapbars;
               
        Bar_Axes_h = axes('Position', BarPosition);
        ColorBarVector = [linspace(color_scale(1),color_scale(2),128)]';      % makes a vector of 128 points
        imagesc([0 1],ColorBarVector,[ColorBarVector ColorBarVector]);                              % scales color in the image to full range
        set(Bar_Axes_h,'XTick',[],'fontsize',12, 'fontweight','bold','Ydir','normal');
        
        %-----  Second Plot
   
        Plot2_Axes_h = axes('position', Position_2);
        m_proj('Equidistant Cylindrical','lon',PlotAxesRange([1 2]),'lat',PlotAxesRange([3 4])); 
        %m_gshhs_h('patch',[.7 .7 .7]);
        m_usercoast('ZK_Coastline', 'patch', LandColor)        % Read from file     
        hold on;
        sub=find(finite(param));
        xplot = lon(sub);
        yplot = lat(sub);
        h_mp = m_plot(xplot, yplot,'b.','markersize',4);
        hold on      
        set(Plot2_Axes_h,'Xtick',[],'Ytick',[]);
        axis(SaveAxis);                                             % Use the same scale as in Plot1
        g=get(Plot2_Axes_h, 'position'); 
        Plot2_Axes_Over_h = axes('position',g);      
        axis(PlotAxesRange); 
        set(Plot2_Axes_Over_h, 'color', 'none', 'visible', 'off'); 
        gry=[.7 .7 .7]; 
        yllocright=0; 
        lblaxes=1; 
        lcol='k'; 
        scfac=.7;       
        mapbars; 
        
        %---- Test here whether jpeg writes properly
        filename=['ZK_TEST.jpg;'];       
        eval(['print -djpeg ',filename]);    
              
        %------------------------------
        AX=[min(lon) max(lon) min(lat) max(lat)];
        
      
         
        %----------------------------------------------------------    
     
         
        %pcolor([0 1],cb,[cb cb]); shading flat;
        %set(gca,'XTick',[],'fontsize',12, 'fontweight','bold');
        
        %if count==1
        %    set(gca,'YTick',CBT,'YTickLabel',CBTL);
        %end

        if count==2
            set(gca,'YTick', ColorBarTick2, 'YTickLabel', ColorBarLabel2);
        end
        
        frame;
        if season_chk == 5,
         t=text(.5,.96,['Depth: ',depthn,'m, Average ',sname ,' ', char(pname(count)), ' Concentrations: ', int2str(min(time)), ':', int2str(max(time))]);
        end
        if (season_chk==1 |season_chk==2 |season_chk==3 |season_chk== 4),
         t=text(.5,.96,['Depth: ',depthn,'m, Southern ',sname,', Average ', char(pname(count)), ' Concentrations: ', int2str(min(time)), ':', int2str(max(time))]);
        end
        set(t,'horizontalalignment','center','fontsize',12,'fontweight','bold');
        t=text(P2(1)+.215,P2(2)+.76,'Station Locations'); set(t,'fontsize',12,'fontweight','bold');
        filename=['peru_contours_', paramss(count,:),'_',sname,'_',depthn, '.jpg;'];
        
        eval(['print -djpeg ',filename]);    
     disp ('end')   
    end %for count = 1:3
  end % for seasons