Vous êtes ici : Accueil / Alexandre Ganachaud / xbt / cmpxbttsg / xbttsg_filter.m

xbttsg_filter.m

Par Webmaster Legos Dernière modification 03/07/2012 18:30

Objective-C source code icon xbttsg_filter.m — Objective-C source code, 3 kB (3621 bytes)

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

Actions sur le document

logo cnes logo IRD Logo université de Toulouse Logo université Paul Sabatier Logo CNRS
Logo bibliothèque OBS Logo Observatoire Midi Pyrénées