xbttsg_filter.m

Contenu du fichier

% NOW ELIMINATE BAD POINTS BASED ON XBT STD DEVIATION BY LAT BANDS
%SEPARATE CRUISES
ca
load Data2/thxbt_corrected

%LATITUDE BANDS FOR STD DEVIATION
%glatband=[-90 -50 -35 -25 -15 -6 3 15 25 35 50 90];
glatband=[-90:5:90];
nstd=2; % mark points out of nstd std deviations;
tthresh=1; %then keep those points if within tthresh deg. of tsg

xbtsurf=[];xlat=[];
%Statistics over all xbt any ship by latitude bands
for iship =1: length(thxbt)
  nxbt=length(thxbt(iship).xbt);
  for ix=1:nxbt
    %All average surface XBT measurements in xbtsurf
    xbtsurf=[xbtsurf mean(thxbt(iship).xbt(ix).xsst)];
    xlat=[xlat thxbt(iship).xbt(ix).xlat];
  end
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculates statistics for all ships by latitude bands
for iband=1:length(glatband)-1
  latmin=glatband(iband);
  latmax=glatband(iband+1);
  gitirs=find(xlat>latmin & xlat <=latmax);
  if ~isempty(gitirs) & length(gitirs > 4)
    stdxbt(iband)=std(xbtsurf(gitirs));
    mnxbt(iband)=mean(xbtsurf(gitirs));
  else
    stdxbt(iband)=NaN;
    mnxbt(iband)=NaN;
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%For each ship Find points out of 3 std deviations
%But keep point if temperature close to TSG measurement (0.8deg.)
%Then correct position of XBT using TSG
ip=0;figure(1);clf
badshots=[]; %ID of bad shots
notsobadshots=[]; %ID of shots with position corrections
for iship =1: length(thxbt)
  clear xbtsurf xlat
  ip=ip+1;
  thg=thxbt(iship).thg;
  xbt=thxbt(iship).xbt;
  gibad=[];
  nxbt=length(xbt);
  tlat=cat(1,thg.thlat);
  tlon=cat(1,thg.thlon);
  xlat=cat(1,xbt.xlat);
  xlon=cat(1,xbt.xlon);
  for ix=1:nxbt
    xbtsurf(ix)=mean(xbt(ix).xsst);
  end
  %plot(xlat,xbtsurf,'x')
  % mark points out of nstd std deviations;
  for iband=1:length(glatband)-1
    latmin=glatband(iband);
    latmax=glatband(iband+1);
    gitirs=find(xlat>latmin & xlat <=latmax);
    gibadl=find(xbtsurf(gitirs)<(mnxbt(iband)-nstd*stdxbt(iband)) | ...
                xbtsurf(gitirs)>(mnxbt(iband)+nstd*stdxbt(iband)));
    gibad=[gibad gitirs(gibadl)'];
    %mark points without statistics
    if isnan(stdxbt(iband))
      gibad=[gibad gitirs];
    end
    %TSG
    %HERE SHOULD KEEP BAD POINTS WITH THG OK
  end %on iband
  %Check if point close to TSG measurement
  %If not remove point
  gisnotsobad=find(abs(cat(1,thg(gibad).thsst)-xbtsurf(gibad)')<0.8);
  g2chg=gibad(gisnotsobad);

  %RECORD ID OF BAD SHOTS
  notsobadshots=[ notsobadshots xbt(g2chg).xtirid];
  badshots=[badshots xbt(gibad).xtirid];
  
  gibad(gisnotsobad)=[];
  thxbt(iship).thg(gibad)=[];
  thxbt(iship).xbt(gibad)=[];

  %Plot what remains in blue, and mark bad points
  clf;figure(gcf)
  plot(xlat,xbtsurf,'+')
  hold on;
  plot(xlat(gibad),xbtsurf(gibad),'ro')
  title(['Corrections on ' thxbt(iship).shipid])
  ppause
end %iship



%PLOT WHAT REMAINS WITH "BAD" XBT REMOVED
ip=0;figure(1);clf
for iship =1: length(thxbt)
  ip=ip+1;
  subplot(2,2,ip)
  thg=thxbt(iship).thg;
  xbt=thxbt(iship).xbt;
  nxbt=length(xbt);
  for ix=1:nxbt
    plot(thg(ix).thsst-xbt(ix).xsst,xbt(ix).xdep,'.');hold on
  end %ix
  xlabel('THG - XBT');
  set(gca,'ydir','rev');ylabel('depth(m)')
  title(sprintf('%s %s-%s',thxbt(iship).shipid,...
    datestr(thg(1).thdate,12),datestr(thg(nxbt).thdate,12)))
  if ip==4
    ip=0;
    drawnow;
    figure(gcf+1);clf
  end
end %on iship

save Data2/hxbt_filtered.mat thxbt badshots notsobadshots
%badshots = away 3 std and away TSG
%notsobadshots =  away 3 std and with .8 deg TSG