% plot_sat_all.m
%Version Dec,2005
% ***this routine created in matlab 6.5 (windows version) - mac users with 7.1 have problems
%     with colorbar****
%
% Near-real time SST data are from microwave sensors. The units are degrees
%  C and are on 1/4 degree grid. 
%
% Near-real time SSH data consist of anomaly data from altimeters plus a
%  mean developed by Kathryn Kelly and Shenfu Dong at the University of 
%  Washington from Topex/Poseidon plus Hydrobase data, please see
%  (http://ultrasat.apl.washington.edu/kkelly/natl). Anomaly data come from
%  the previsou ten days of measurements from the Jason and Geosat Follow On 
%  (GFO) altimeters, courtesy of Colorado Center for Astrodynamics Research, 
%  (CCAR). The units are meters and the data are on a 1/4 degree grid.
%
% Near-real time surface wind vectors come from the QuikSCAT scatterometer.
%  The wind file has several variables. Wind speed (wspd), wind components 
%  (u and v) and minutes (min) in the day (UTC) when the data were taken. 
%  Each variable is separated into two maps (to preclude overlap). One map 
%  consists of the ascending orbits (SE to NW) where the variable name 
%  ends in "a", ie wspda and the other contains the data from the descending 
%  orbits (NE to SW) where the variable name ends in "d"ie, wspdd. The time 
%  between successive orbits is about 100 minutes. The time between ascending 
%  orbits and descending orbits is roughly half a day. There are gaps between 
%  the orbits which show up as NaNs in the .mat file. If the rain flag is 
%  postive and the wind speed is less than 10m/s, we replace the
%  data with a NaN. The winds are in m/s and are gridded to 1/4 degree.
%  S.Dickinson Sept 2005
%
%23 May 2005 - jadb
% This version is a composite of Suzannes'  load_nrt_sst/ssh/
% It assumes standard file naming convention [var]yyyymmdd.mat,
%   where var can be 'sst' 'ssh' or 'win' without the quotes.
% It loads gridded data from same directory and plots
%  using coast.mat from \toolbox\map\mapdisp to plot coastline
%
%28 Sept 2005 - add x/y extraction option. Allows user to
% locate waypoints on an image (usually sst). first zoom and/or
% hit return. then put in waypoints (queried if okay after each one)
% then hit return to terminate. lon lat position are listed. 
% add station position between way points.- jadb
%
%30 sept 2005 - change variable names in wind file
%   wspd_am=wspda...etc - jadb
%
%Nov 2005 - add terry's geovel plot
%           sta positions to deg & min
%Dec 2005 - add o/p file for storing positions, and to give to bridge. -
%jadb
%
clear
close('all','hidden')
%
disp(['plot_sat_all   ver dec05']);
disp(' ');disp(' ');
%
statdd=[];
statdm=[];
dowhich=[];
getpos='n';
docont='n';
dotrack='n';
domoor='n';
dobar='n';
hflag='n';
isasis='n';
defvec=[17.5 18.5];
defbounds=[-80 -50 30 50];
posfile='knorr.pos'
%
%
popt=['  1 = sst            ';...
      '  2 = sst w/ velocity';...
      '  3 = win_am         ';...
      '  4 = win_pm         ';...
      '  5 = ssh            '];

disp([' plotting options:  ']);
disp(popt);
dopop=input(' Please enter selection:  ');
switch dopop
    case 1
        dtype='sst';
    case 2
        dtype='ssh';  
        dowhat='vec';
    case 3,
        dtype='win';
        dowhich='m';
    case 4
        dtype='win';
        dowhich='e';
    case 5
        dtype='ssh';
        dowhat='con';
 end%switch
%           
ls *.mat
fdate=input('Enter date (ie 20051107):  ','s');
filein=[dtype,fdate,'.mat'];
eval(['load ',filein]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if dtype == 'sst'
  docont=input('Plot sst contours(n=def)?   ','s');
  if lower(docont) == 'y'
      tconturs=input('Enter 2 contour intervals (default=  [17.7 18.5]  ):  '); 
      if isempty(tconturs)
          tconturs=defvec;
      end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bounds=input('Enter bounds: ie [-80 -50 30 50]   ');
if isempty(bounds)
    bounds=defbounds;
end
if max(lon) > 200;lon=lon-360;end

xmin=min(lon);xmax=max(lon);ymin=min(lat);ymax=max(lat);
flat=find(lat>= bounds(3) & lat <= bounds(4));
flon=find(lon>= bounds(1) & lon <= bounds(2));
lat=lat(flat);
lon=lon(flon);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dobar=input('Colorbar? (n=def):  ','s');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
domoor=input('Overplot mooring locations(n=def)?   ','s');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dotrack=input('Overplot longitude latitude tracks(def=n)?   ','s');
if lower(dotrack) == 'y'
    k=0;
    tsym=['*','+','s','o'];
    tline=['-','--','-.',':'];
    tcolr=['k','b','w'];
    ntrk=input('enter number of lon/lat files:  ');
    disp([' expecting input file(n x 2) (longitude latitude)']);
    nknt=0;
    for nn=1:ntrk
      nknt=nknt+1;
     disp([' ']);
     tfile=input(' enter input trackfile name:  ','s');
     k=strfind(tfile,'asis');
     if k == 1
%         eval(['load ',tfile]);
%         alon=asis.lon';
%         alat=asis.lat';
         eval(['asis=load(','tfile',')']);
         alon=asis.lon;
         alat=asis.lat;
         trdata(nknt:nknt+length(alon)-1,1)=alon;
         trdata(nknt:nknt+length(alat)-1,2)=alat;
     else
       temp=load(tfile);
       trdata(nknt:nknt+length(temp)-1,1)=temp(:,1);
       trdata(nknt:nknt+length(temp)-1,2)=temp(:,2);
     end
     trtype(nn)=input('symbol or line (s/l) :','s');
     nknt=length(trdata)+1;
     trdata(nknt,1:2)=[999 999];
    end%nn
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
getpos=input('Extract XY cordinates(n=def)?  ','s');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fdot=find(filein == '.');
y=str2num(filein(4:7));
m=str2num(filein(8:9));
d=str2num(filein(10:11));
mon=['Jan';'Feb';'Mar';'Apr';'May';'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec'];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% imagesc plots pixel center...
%
clf
figure(1)
switch dtype
    case 'sst'  
       if ~exist('sst','var')
           error(' gridded sst variable does not exist');
       end
       data=sst(flat,flon);
       tstr=sprintf('SST from %s %d, %d',mon(m,:),d,y);
       imagesc(lon,lat,data);
       hold on
       if lower(docont) == 'y'
          [a,b]=contour(lon,lat,data,tconturs,'k');
          clabel(a,b);
          set(b,'linewidth',2)
       end
       colorbar;axis xy;
   case 'ssh'
       if ~exist('ssh','var')
           error(' gridded ssh variable does not exist');
       end
       data=ssh(flat,flon);
       latg=lat;
       long=lon;
       if dowhat == 'con'
         tstr=sprintf('SSH from %s %d, %d',mon(m,:),d,y);
         imagesc(lon,lat,data);
         colorbar;axis xy;
       elseif dowhat == 'vec'
 %new file with different bounds
         filein2=['sst',filein(4:end)];
         if ~exist(filein2,'file')
              error(' sst file for this date does not exist');
         end
           eval(['load ',filein2])
           if max(lon) > 200;lon=lon-360;end
           flat=find(lat>= bounds(3) & lat <= bounds(4));
           flon=find(lon>= bounds(1) & lon <= bounds(2));
           lat=lat(flat);
           lon=lon(flon);
           data2=sst(flat,flon);
         tstr=sprintf('SST & surface velocity for %s %d, %d',mon(m,:),d,y);
  % this from terry's geovel
         [dx0 dy0]=lalon(mean(lat));dx0=1000*dx0;dy0=1000*dy0; %dist. in meters
         yg=dy0*(latg-mean(latg));
         xg=dx0*(long-mean(long));
         deltax=mean(diff(xg));
         deltay=mean(diff(yg));
         [vg,ug] = GRADIENT(data,deltax,deltay); %unscaled velocities
            ug=-ug;
            [iy ix]=size(ug);
         for i=1:iy
             u(i,:)=(sw_g(latg(i),0))/sw_f(latg(i))*ug(i,:);
             v(i,:)=(sw_g(latg(i),0))/sw_f(latg(i))*vg(i,:);
         end
         imagesc(lon,lat,data2,[2 27])
           hold on
         [xl,yl]=meshgrid(long,latg);
          quiver(xl,yl,u(1:1:end,1:1:end),v(1:1:end,1:1:end),3);
         axis xy
         if lower(dobar) == 'y';colorbar;end
     end %if
   case 'win'
        [mlon,mlat]=meshgrid(lon,lat);
       if dowhich == 'm'
           data=wspd_am(flat,flon);
           u=u_am(flat,flon);
           v=v_am(flat,flon);
           tstr=sprintf('Wind speed from %s %d, %d 0900-1000 UTC',mon(m,:),d,y);
       elseif dowhich == 'e'
           data=wspd_pm(flat,flon);
           u=u_pm(flat,flon);
           v=v_pm(flat,flon);
           tstr=sprintf('Wind speed from %s %d, %d 2200-2300 UTC',mon(m,:),d,y);
       end
       imagesc(lon,lat,data,[0 25]);
       colorbar;
       axis xy;
       hold on      
       qq=quiver(mlon,mlat,u,v);
end%switch

       title(tstr,'fontsize',12);
       hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if lower(domoor) == 'y'
   mpmoor=[-56.0 38.0;-60.1 36.2];
   surfmoor=[-64.713 38.339];
   hh=plot(mpmoor(:,1),mpmoor(:,2),'kh');
   set(hh,'color',[252/255 252/255 252/255],'markersize',8,'linewidth',2);
   h2=plot(surfmoor(1),surfmoor(2),'kd');
   set(h2,'color',[252/255 252/255 252/255],'markersize',8,'linewidth',2);
%line W moorings, gusto is furthest one out
   disp(['line W moorings in ''wmoor'' ']);
   wmoor=[-69.7169 39.6013;-69.4450 39.2167;-69.1883 38.844;...
           -68.9026 38.4244;-68.6634 38.0733;-68.2632 37.5191];
   [w1,w2]=size(wmoor);
   for w=1:w1
     ww(w)=text(wmoor(w,1),wmoor(w,2),'w','color','k','fontweight','bold','fontsize',15);
     set(ww(w),'color',[252/255 252/255 252/255],'horizontalalignment','center','VerticalAlignment','baseline');
   end
end
%%%%%%%%%%%%%%%%%%%
% coast.mat contains long and lat
load coast
coastlines(1,:)=long';
coastlines(2,:)=lat';
hold on
cc=plot(coastlines(1,:),coastlines(2,:),'w');
set(cc,'Erasemode','none','linewidth',1);
%
%overide data bounds
axis(bounds)

%%%%%%%%%%%%%%%%%%%
if lower(dotrack) == 'y'
    sknt=0;
    lknt=0;
    strt=1;
   for n=1:ntrk
       ft=find(trdata(:,1) == 999);
       tdata=trdata(strt:ft(n)-1,:);
       if trtype(n) == 's'
          sknt=sknt+1;
          ht(n)=plot(tdata(:,1),tdata(:,2),[tsym(sknt) 'k'],'markersize',8');
       elseif trtype(n) == 'l'
          lknt=lknt+1;
          ht(n)=plot(tdata(:,1),tdata(:,2),[tline(lknt) 'k'],'linewidth',2);
       end%trtype
      strt=ft(n)+1;
   end%for
end%if
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if lower(getpos) == 'y'
    %%%%%o/p position file (stp=1 if new posfile, w/append stp=2)
     n=1;
     stp=1;
   if exist(posfile,'file')
           !copy atlantis13.pos atlantis13.bck
     posdata=load(posfile);
     pp=plot(posdata(1:end-1,1),posdata(1:end-1,2),'*k');

     laststa=posdata(end,3);
     x(1)=posdata(end,1);
     y(1)=posdata(end,2);
     tp=plot(x(n),y(n),'x');
     set(tp,'color','k','markersize',8,'linewidth',2);
     n=2;
     stp=2; 
   else
     posdata=[];
     laststa=0;
     hflag='y';
   end%exist
%
    disp([' ']);
    disp(['IMPORTANT:  Leave mouse in figure!!']);
    disp(['    if you need to zoom...do it now, then hit return!!']);
    pause
    
    disp(['Click mouse over way positions (return to exit)']);
%%%%%%%%%establish waypoints
  while n < 99
      [x1,y1]=ginput(1);
    if isempty(x1);break;end
    x(n)=x1;y(n)=y1;
    tdot(n)=text(x(n),y(n),'.','color','k','fontweight','bold','fontsize',25);
    set(tdot(n),'horizontalalignment','center','VerticalAlignment','baseline');
    tcirc(n)=text(x(n),y(n),'O','color','y','fontsize',20);
    set(tcirc(n),'horizontalalignment','center','VerticalAlignment','middle');

    isokay=input(' Is this OKAY(y/n)?  ','s');
      if lower(isokay) == 'y'
         n=n+1;
      else
        delete(tdot(n));
        delete(tcirc(n));
      end
  end%while
  disp([' ']);
  for k=1:length(x)
    disp([' waypoint ',num2str(k),': ',num2str(y(k)),' ',num2str(x(k))]);
  end 
%%%%%%%%%%%%%%%%%% stations between way points
 disp([' ']);
 docalc=input(' Calculate station positions between end points (y/n)?:  ','s');
 disp([' ']);
 if lower(docalc) == 'y'
   for k=1:length(x)-1
       [xdist,xang]=sw_dist(y(k:k+1),x(k:k+1),'km');
    disp(['Between waypoints ',num2str(k),' and ',num2str(k+1),'   ',num2str(round(xdist)),' km']);
    nums=input('   enter number of stations: ');
    if nums > 0
      xdiff=(x(k+1)-x(k))/(nums+1);
      for nn=1:nums
         xx(nn)=x(k) + nn*(xdiff);
         yy(nn)=intrplat(x(k:k+1),y(k:k+1),xx(nn));
      end
      statpos(:,1)=[y(k);yy';y(k+1)];
      statpos(:,2)=[x(k);xx';x(k+1)];
  %convert to deg & min
   [nr,nc]=size(statpos);
  % if more that one waypoint, need to avoid duplication (stp=start)
  statdd=[statdd;statpos(stp:end,[2 1])];
  % try to add in sta numbers
    for r=1:nr
       dpos(r,1)=floor(statpos(r,1));
       dpos(r,2)=(statpos(r,1)-dpos(r,1))*60;
       dpos(r,3)=(floor(abs(statpos(r,2))));
       dpos(r,4)=(abs(statpos(r,2))-dpos(r,3))*60;
       dpos(r,3)=dpos(r,3)*-1;
       outpos(r,:)=sprintf('%5d %7.4f %6d %7.4f',dpos(r,:));
   end%end r
   statdm=[statdm;outpos(stp:end,:)];
   disp(' ');
   disp(outpos);
    disp(' ');
   clear xx clear yy clear xdiff statpos outpos dpos
   stp=2;
  end%nums>0
 end%end k
%need to put sta num in statdm
    [ndr,ndc]=size(statdm);
    statdm2=str2num(statdm);
    for dr=1:ndr
        laststa=laststa+1;
        statdm3(dr,:)=[laststa statdm2(dr,:)];
    end
  
  newdata=[statdd statdm3];
  [nst,npo]=size(newdata);
  for c=1:nst
      cposdata(c,:)=sprintf('%8d %5d %7.4f %6d %7.4f',newdata(c,3:7));
  end
  fp3=fopen(posfile,'a'); 
  if hflag == 'y'
    fprintf(fp3,'%%atlantis 13 \n');
    fprintf(fp3,'%%  lon      lat       sta#  latdeg latmin  londeg lonmin\n');
  end
  for ss=1:nst
      fprintf(fp3,'%7.4f  %7.4f %37s\n',newdata(ss,1:2),cposdata(ss,:));
  end
  fclose(fp3);
   disp([' positions appended to file: ',posfile]);
  end%docalc
end%getpos

%