xbttsg_filter.m
by
Webmaster Legos
—
last modified
Jul 03, 2012 06:30 PM
xbttsg_filter.m
—
Objective-C source code,
3 kB (3,621 bytes)
File contents
% 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