
% MATLAB routine to convert various satellite images to MATLAB contours

% Note that many of the prompts are commented out (i.e. hardcoded)
% Note that one may need to access to special mfiles in the "atsee" directory


%clf


% yr=input('What year 2-digits?');
yr=99;
disp(' ');
disp('Harcoded Assumptions listed in this window.');
disp(' ');
disp(['Assumes year = ',num2str(yr)]);


%hmi=input('How many images per page?');
hmi=2;
disp(' ');


% batch mode file read
fid1=fopen('pthname.dat','r');
%hmf=input('How many files in your "pthname.dat" list? ');
hmf=hmi;
disp(['Assumes you are processing the latest images file listed in "pthname.dat"']);

% increment of resolution (ie every "inc" pixels used)
%inc=input('Use every Nth point .. enter N: ');
inc=1;

mmm=['Jan';'Feb';'Mar';'Apr';'May';'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec'];
wtype=['revd';'rmp ';'nav ';'pvu ';'byt ';'merg';'ns7 ';'g*7 ';'gbs ';'dim2'];

%%%%%%%%%%%%%%%%%%  BIG LOOP THROUGH IMAGES START HERE %%%%%%%%%%%%%%%%%

for iii=1:hmf



% only make "subplots" sometimes when the # of images >1
  if hmi~=1;
%    subplot(1,2,iii);
  else;
%    clf;
  end;

  % wt=menu('what file type?','revd','rmp','nav','pvu','byt','merg','ns7','g*7','gbs');
  if iii==1;
    wt=9;
    % read next filename from "pthname.dat" file
    [fn]=fscanf(fid1,'%s',1);
    fn=['/data5/jmanning/images/1999/',fn];
    disp(' ');
    disp(['Processing ',fn,'...']);
  elseif iii==2;
    wt=10;
    fn(length(fn)-2:length(fn))='dim';
    fn=[fn,'2'];
  end;
  disp(' ');
  disp(['Assumes you are processing ',wtype(wt,:),' type images.']);

  %  call a subroutine that gets the "counts" from satellite image
  getimage_batch;

  % define lat&lon
  % old options uses entire input field
  if wt==1|wt==6|wt==9;
    lon=-70.632:.01305*inc:-63.955;
    lat=38.994:.00977*inc:43.996;
    mnlo=-70.632;mxlo=-63.955;
    mnla=38.994;mxla=43.996;
    
  elseif wt==2|wt==3;
    mxla=45.42900; mnla=39.02500; mxlo=-63.50400; mnlo=-72.16400;
    ila=(mxla-mnla)/510*inc; ilo=(mxlo-mnlo)/510*inc;
    lon=mnlo:ilo:mxlo;
    lat=mnla:ila:mxla;
    % correction to center pixel?? I.E. MANUAL NAVIGATION OF IMAGE!!!
    lon=lon+ilo/2;
    lat=lat+ila;

  elseif wt==4;
    %    mxla=53.778; mnla= 28.240; mxlo=-54.200; mnlo=-85.366;
    %    mxla=60.575; mnla= 9.499; mxlo=-33.89; mnlo=-96.22;
    ila=(mxla-mnla)/511*inc; ilo=(mxlo-mnlo)/511*inc;
    lon=mnlo:ilo:mxlo;
    lat=mnla:ila:mxla;
  elseif wt==5;
    mxla=42.5; mnla= 39.0; mxlo=-68.5; mnlo=-72.5;
    ila=(mxla-mnla)/684*inc*(-1); ilo=(mxlo-mnlo)/791*inc;
    lon=mnlo:ilo:mxlo;
    lat=mxla:ila:mnla;
  elseif wt==7;
    mxla=45.82; mnla= 29.930; mxlo=-58.81; mnlo=-79.078;
    %   ila=(mxla-mnla)/511*inc*(-1); ilo=(mxlo-mnlo)/511*inc;
    ila=(mxla-mnla)/511*inc*(1); ilo=(mxlo-mnlo)/511*inc;
    lon=mnlo:ilo:mxlo;
    %    lat=mxla:ila:mnla;
    lat=mnla:ila:mxla;
  elseif wt==8;
    mxla=43.71875; mnla=38.640625; mxlo=-65.195312; mnlo=-71.95312;
    %    ila=(mxla-mnla)/511*inc*(-1); ilo=(mxlo-mnlo)/511*inc;
    ila=(mxla-mnla)/511*inc*(1); ilo=(mxlo-mnlo)/511*inc;
    lon=mnlo:ilo:mxlo;
    %    lat=mxla:ila:mnla;
    lat=mnla:ila:mxla;
  end;

  % For the case of frontal locations, skip over next several loops
  if wt~=10;

    % convert to temperature and rotate to frame
    if wt~=5&wt<7|wt==9;
      sst=rot90(0.125*counts(1:inc:512,1:inc:512));
      cloud=4*0.125;
    elseif wt==5;
      sst=counts/10;
      cloud=0;
    elseif wt>=7&wt~=9;
      sst=rot90(counts(1:inc:512,1:inc:512));
      sst=sst/10;
      cloud=0;
    end;
    mingood=3;
    maxgood=20;
    disp(' ');
    disp(['Assumes you are interested in the range: ',num2str(mingood),...
          ' - ',num2str(maxgood),' degrees C']);
    [nn,mm]=size(sst);

    %   get rid of cloud values
    for im=1:mm;
      for in=1:nn;
        if sst(in,im)>=cloud-.05&sst(in,im)<=cloud+.05;
          sst(in,im)=NaN;
        end;
        if sst(in,im)<mingood|sst(in,im)>maxgood;
          sst(in,im)=NaN;
        end;
      end;
    end;

    % this option reduces the field prior to contouring so that
    % "contourfill" and "econtour" can be used
    %        wo=menu('which field?',...
    %                '1 just georges',...
    %                '2 georges include GS-rings (WCRII)',...
    %                '3 GS-rings & gulfstream',...
    %                '4 entire',...
    %                '5 Gulf of Maine',... 
    %                 '6 Coastal Mixing & Optics area',...
    %                '7 just southern flank',...
    %                '8 just GS-rings');
    wo=7;

    if wo==1;
      id1=find(lat>=39.9&lat<=43.1);
      id2=find(lon>=-69.6&lon<=-65.4);
    elseif wo==2;
      id1=find(lat>=39.&lat<=45);
      id2=find(lon>=-71.5&lon<=-65.0);
    elseif wo==3;
      id1=find(lat>=36.0&lat<=43);
      id2=find(lon>=-71&lon<=-64.);
      disp('Hardcoded to just GS ring area [-71 -63 37 43].');
    elseif wo==4;
      id1=find(lat>=mnla&lat<=mxla);
      id2=find(lon>=mnlo&lon<=mxlo);
    elseif wo==5;
      id1=find(lat>=40&lat<=45.5);
      id2=find(lon>=-71&lon<=-65);
    elseif wo==6;
      id1=find(lat>=39&lat<=42.5);
      id2=find(lon>=-72.5&lon<=-68.5);
    elseif wo==7;
      id1=find(lat>=40.&lat<=42.);
      id2=find(lon>=-68.&lon<=-66.5);
    elseif wo==8;
      id1=find(lat>=39.0&lat<=42.5);
      id2=find(lon>=-70&lon<=-65.);
      disp('Hardcoded to just GS ring area [-70 -65 39 42].');
    end

    sst=sst(id1,id2);

    %  try forcing min&max one data
    %  sst(1,1)=4;
    %  sst(1,2)=24;
    hold on;

    % report range of sst
    [rg,mint,maxt]=range(sst);

    % Now the data is loaded .. start the plotting phase
    x=lon(id2);y=lat(id1);z=sst;

    % heres my attempt to prep data for contouring
    % grdit=menu('Do you want to grid data?','y','n');
    grdit=2;
    disp(' ');
    disp('Assumes no gridding of data.');
    if grdit==1;
         [nr,nc]=size(z);
         disp('WARNING: May need to review the grid setup');
    %      for isat=1:nr;
    %        xsat=[xsat x];
    %      end;
    %      for isat=1:nc;
    %        ysat=[ysat y];
    %      end;
         [X,Y]=meshgrid(x,y);
         xsat=reshape(X,1,nr*nc);
         ysat=reshape(Y,1,nr*nc);   
         zsat=reshape(z,1,nr*nc);
         id=find(isnan(zsat)~=1);
         x=xsat(id)';y=ysat(id)';z=zsat(id)';
    end;

    hsec_batch


    if wo==1;
      axis([-69 -66 40 43]);
    elseif wo==2;
      axis([-68.5 -67. 39.5 41.5]);
    elseif wo==3;
      axis([-69.5 -64.5 37.5 42.5]);
    elseif wo==4;
      axis([mnlo mxlo mnla mxla]);
    elseif wo==5;
      axis([-71 -65.5 41 45]);
    elseif wo==6;
      axis([-72.5 -68.5 39 42.5]);
    elseif wo==7;
      axis([-68. -66.5 40.5 41.5])
    elseif wo==8;
      axis([-69.5 -65.5 39.75 41.25])
    end;
    set(gca,'TickDir','out');

  end;
  % processing the fronts now as well


  % change the filename starting at last slash
  ids=find(fn=='/');
  if ids~=[];
    fn2=fn(ids(length(ids))+1:length(fn));
  end;
  if wt>=7;
    lh=2;
  else;
    lh=4;
  end;
  sd=4;
  if fn2(length(fn2)-7:length(fn2)-5)=='RAS';
    cdy=calday(yr,str2num(fn2(length(fn2)-4:length(fn2)-2)));
    month=mmm(cdy(:,2),:);
    title([num2str(cdy(:,3)),' ',month,' 19',num2str(yr),' @',fn2(length(fn2)-1:length(fn2)),'h gmt'],'fontsize',16);
  else;
    yr=str2num(fn2(2:3));
    if yr>75&yr<100;
      cdy=calday(yr,str2num(fn2(sd:sd+2)));
      month=mmm(cdy(2),:);
      title([num2str(cdy(:,3)),' ',month,' 19',num2str(yr),' @',fn2(7:8),'h UTC'],...
      'color','r','fontsize',16);
    else
      title(fn2,'fontsize',12);
    end;
  end;

  if iii==1;
    % bathy contour labels
    text(-69.6,40.55,'60m','rotation',-15,'fontsize',12,'color','k','horizontalalignment','right')
    %  text(-69.6,40.33,'80m','rotation',-25,'fontsize',12,'color','k','horizontalalignment','right')
    text(-69.2,40.1,'100m','rotation',10,'fontsize',12,'color','k','horizontalalignment','right')
    text(-69.2,39.9,'200m','rotation',5,'fontsize',12,'color','k','horizontalalignment','right')
  end;

  % add other postings (drifters, moorings, ctd stations, etc.)
  %apost=menu('Will you be adding posted data?','y','n');
  apost=1;

  if apost==1;
    disp('Harcoded to  post data.');
    %    pwhich=menu('Post What?','1995 drifters','wcr','broadscale stas',...
    %                 '1995 moorings','1997 drifters','1997 Rons',...
    %                 'may95 crutrk','1997 drifters & mocness',...
    %                 '1999 Rons Moorings');
    pwhich=9;
    if pwhich==1;
    %     plot '95 drifter deployments
    %     load /data1/jgofs/sj9507/drift/p1a9507e.mat -ascii
    %     d1=p1a9507e;
    %     plot(d1(:,4),d1(:,3),'k-','linewidth',2)
      sumplt;

    elseif pwhich==2;
      path(path,'/home3/ocn/jmanning/wcr/mfiles');
      overlaysat;

    elseif pwhich==3;
      pltbs
    elseif pwhich==4;
      %     process
      %      plot([-67.6264 -67.5588],[40.956 40.864],'r.','markersize',20);
      %     longterm
      %      plot([-67.317 -67.542],[40.968 41.407],'r.','markersize',20);
      %     canadian
      %      plot(-66.5347,41.7335,'r.','markersize',20);
      %      fvelip(str2num(fn(4:6)));
    elseif pwhich==5;
      sumplt;
    elseif pwhich==6|pwhich==9;
      hold on;
      pltrons(2);
      if pwhich==9;
        text(-67.25,41.17,'Moorings ','rotation',-65,'color','y');
      end;
    elseif pwhich==7;
      load /data5/jmanning/ship/sj9507/sj9507
      plot(-y(:,3),y(:,2),'k');
    elseif pwhich==8;
      overmoc;
      pv=menu('plotvel?','y','n');
      if pv==1;
        plotvel;
      end;
    end;

  end;


  % position and label according to subplot
  if hmi==3;
    if iii==1;
      set(gca,'position',[.1 .1 .28 .8]);  
      %      set(gca,'xticklabels',[' '])
      %      xlabel(' ');
    elseif iii==2;
      set(gca,'position',[.4 .1 .28 .8]);
      %      set(gca,'xticklabels',[' '])
      %      xlabel(' ');
      set(gca,'yticklabels',[' '])
      ylabel(' ');
    elseif iii==3;
      set(gca,'position',[.7 .1 .28 .8]);
      set(gca,'yticklabels',[' '])
      ylabel(' ');
    elseif iii==4;
      set(gca,'position',[.5 .08 .4 .4]);
      set(gca,'yticklabels',[' '])
      ylabel(' ');
    end;
  end;

  if hmi==2&batch~=1;
    if iii==1;
      set(gca,'position',[.1 .1 .425 .8]);  
      %      set(gca,'xticklabels',[' '])
      %      xlabel(' ');
    elseif iii==2;
      set(gca,'position',[.55 .1 .425 .8]);
      %      set(gca,'xticklabels',[' '])
      %      xlabel(' ');
      ylabel(' ');
    end;
  end;

  if exist('fid')~=0;
    fclose(fid);
  end;
  drawnow;
  hold on;

  % do next file/image
end;

fclose(fid1);

%set(gcf,'paperposition',[.5 1 8 9.5]);

if hmi>2&hmi<4;
  orient landscape;
else;
%  orient tall;
end;

eval(['print -dgif8 /data5/jmanning/images/1999/',month,num2str(cdy(3)),'@',fn2(8:11),'Z']);
!/usr/local/bin/elm -s new_sst_imagery jmanning < /home3/ocn/jmanning/gbs99/sat/newimage.txt

