% plot up some tidal and low-frequency stats

figure % plot m2 major axis current amplitude
plot(m2maj_obs(:),m2maj_mod(:),'o',[0 1],[0 1],'k-');
grid;axis square
xlabel('obs');ylabel('mod');title('M2 tide major axis amplitude (m/s)')

figure % plot m2 current phase
plot(m2pha_obs(:),m2pha_mod(:),'o',[-50 50],[-50 50],'k-');
grid;axis square
xlabel('obs');ylabel('mod');title('M2 tide phase');

% plot mean flow from Data and Model
bbox=sscanf(q.bbox,'%g,%g,%g,%g'); % get lon/lat from opensearch request
structure.lon=bbox([1 3]);
structure.lat=bbox([2 4]);
structure.xy_stride=[ 2 2];
%topo = ncgeodataset('c:\rps\cf\fvcom\gom15.nc');
topo = ncgeodataset('http://geoport.whoi.edu/thredds/dodsC/bathy/gom15');
tvar = topo.geovariable('topo');
subs = tvar.geosubset(structure);
load gom_coast.mat
%%
figure % mean flow vectors: DATA
ind=find(subs.data>0);
subs.data(ind)=nan;
pcolorjw(subs.grid.lon,subs.grid.lat,subs.data);
cax=caxis;cax(1)=max(cax(1),-400);caxis(cax);

line(coast(:,1),coast(:,2),'color','black')
hold on;
plot(lon_obs,lat_obs,'ko')
a=arrows(lon_obs,lat_obs,mean1,5,'black')
al=arrows(-71,44.5,0.1,5,'black');
text(-71,44.6,'0.1 m/s')
title('Mean Flow: DATA')
dasp(44)
%%
figure % mean flow vectors: MODEL
pcolorjw(subs.grid.lon,subs.grid.lat,subs.data);
cax=caxis;cax(1)=max(cax(1),-400);caxis(cax);
line(coast(:,1),coast(:,2),'color','black')
hold on;
plot(lon_mod,lat_mod,'ko')
a=arrows(lon_mod,lat_mod,mean2,5,'black')
al=arrows(-71,44.5,0.1,5,'black');
text(-71,44.6,'0.1 m/s')
title('Mean Flow: MODEL')
dasp(44)

figure % plot low-frequency complex correlation between MODEL/DATA as colored dots
pcolorjw(subs.grid.lon,subs.grid.lat,subs.data);
cax=caxis;cax(1)=max(cax(1),-400);caxis(cax);
line(coast(:,1),coast(:,2),'color','black')
hold on;
plot(lon_mod,lat_mod,'ko')
cdot(lon_obs,lat_obs,abs(corr),jet,20,0,[0 1])
title('DATA/MODEL low-frequency correlation: blue=0, green=0.5, red=1.0')
dasp(44)


% interpolate water depth from 3 sec bathymetry grid
%url='c:/rps/cf/fvcom/gom03_v07.nc';
url='http://geoport.whoi.edu/thredds/dodsC/bathy/gom03_v03';
nc=ncgeodataset(url);
lonb=nc{'lon'}(:);
latb=nc{'lat'}(:);
for imoor=1:length(lon_obs);
  i=near(lonb,lon_obs(imoor));
  j=near(latb,lat_obs(imoor));
  wdepth(imoor)=nc{'topo'}(j,i);
end

for imoor=1:length(lon_obs);
  nc=ncgeodataset(dap{imoor});
  wd=nc.attribute('WATER_DEPTH');
  if isempty(wd)
    wdepth2(imoor)=nan;
  else
    wdepth2(imoor)=wd;
  end
end
% pull off names of moorings from DAP urls
for imoor=1:length(lon_obs);
  parts=regexp(dap{imoor},'/','split');
  parts=regexp(parts(end),'\.','split');
  iname{imoor}=parts{1}(1);
end
%% plot low frequency stick plots and create KML

% read header, trailer and point template for KML construction
fid=fopen('header.kml');
str_header=fscanf(fid,'%c');
fclose(fid)

fid=fopen('trailer.kml'); % read trailer of output file to append later
str_trailer=fscanf(fid,'%c');
fclose(fid)

fid = fopen('point_template.kml');  % read KML template for each mooring
f = fscanf(fid, '%c');
fclose(fid);


figure;
set(gcf,'pos',[20 50 1000 700])
set(gcf,'color','white');

fout=fopen('month.kml','w'); % open output KML file
fprintf(fout,'%s\n',str_header); % write header to output file

for imoor=1:length(lon_obs);
%for imoor=8:19;
  %imoor=1;
  smax=max(max(abs(wmodlp{imoor}(:)),max(abs(wobslp{imoor}(:)))));
  vmax=max(max(imag(wmodlp{imoor}(:)),max(imag(wobslp{imoor}(:)))));
  vmin=min(min(imag(wmodlp{imoor}(:)),min(imag(wobslp{imoor}(:)))));
  outf=sprintf('plot_%3.3d',imoor);
  h=timeplt(jdmat2jdrps(jdlp{imoor}),...
    [wobslp{imoor} wmodlp{imoor} abs(wobslp{imoor}) abs(wmodlp{imoor})],...
    [2 1 3 3],...
    [vmin vmax;vmin vmax;0 smax]);
  stacklbl(h(1),'MODEL','m/s');
  stacklbl(h(2),'DATA','m/s');
  stacklbl(h(3),'SPEED','m/s');legend('DATA','MODEL');
  titl=sprintf('Low-Frequency (<33hr) Current Comparison at %s: lon=%7.3f,lat=%7.3f,depth=%6.1f,water depth=%6.1f\n',...
    char(iname{imoor}),lon_obs(imoor),lat_obs(imoor),dep_obs(imoor),wdepth(imoor));
  title(titl,'interpreter','none');
  fixpaper2
  eval(['print -dpng ./frames/' outf]);
  % write kml
  
  
  f2=f;
  f2=strrep(f2,'@NAME@',sprintf('%s',char(iname{imoor})));
  local=0
  if local
    f2=strrep(f2,'@URL@',...
    sprintf('file:///C:/RPS/cf/fvcom/frames/%s.png',outf));
  else
    f2=strrep(f2,'@URL@',...
    sprintf('http://coast-enviro.er.usgs.gov/models/share/frames/%s.png',outf));
  end
  f2=strrep(f2,'@LAT@',sprintf('%12.6f',lat_obs(imoor)));
  f2=strrep(f2,'@LON@',sprintf('%12.6f',lon_obs(imoor)));
  f2=strrep(f2,'@START@',sprintf('%s',datestr(jd{imoor}(1),'yyyy-mm-ddTHH:MM:SSZ')));
  f2=strrep(f2,'@STOP@',sprintf('%s',datestr(jd{imoor}(end),'yyyy-mm-ddTHH:MM:SSZ')));
  f2=strrep(f2,'@WDEPTH@',sprintf('%6.2f',wdepth(imoor)));
  f2=strrep(f2,'@IDEPTH@',sprintf('%6.2f',dep_obs(imoor)));
  fprintf(fout,'%s\n',f2)
  
end
fprintf(fout,'%s',str_trailer);