VMT_RepBadGPS

PURPOSE ^

Replaces bad GPS data with bottom track data.

SYNOPSIS ^

function A = VMT_RepBadGPS(z,A)

DESCRIPTION ^

 Replaces bad GPS data with bottom track data.

 (adapted from code by J. Czuba)

 Known bugs--Will result in bad values if GPS and Bottom track are both
 bad.  FIXED-9-8-09 with interpolation of missing bottom track from good
 data (may have issues with lots of missing data).

 Added capability to detect fly-away GPS points. Works in most cases, but
 may still miss a few fly-aways. Identified fly-aways are replaced with
 interpolation of bottom track (FLE 12-12-12)

 Rearranged the identification of bad values to use logical indexing,
 which is faster, and easier to trace (FLE 12-12-12)

 P.R. Jackson, USGS, 12-9-08 
 Last modified: F.L. Engel, USGS 12/12/2012

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function A = VMT_RepBadGPS(z,A)
0002 % Replaces bad GPS data with bottom track data.
0003 %
0004 % (adapted from code by J. Czuba)
0005 %
0006 % Known bugs--Will result in bad values if GPS and Bottom track are both
0007 % bad.  FIXED-9-8-09 with interpolation of missing bottom track from good
0008 % data (may have issues with lots of missing data).
0009 %
0010 % Added capability to detect fly-away GPS points. Works in most cases, but
0011 % may still miss a few fly-aways. Identified fly-aways are replaced with
0012 % interpolation of bottom track (FLE 12-12-12)
0013 %
0014 % Rearranged the identification of bad values to use logical indexing,
0015 % which is faster, and easier to trace (FLE 12-12-12)
0016 %
0017 % P.R. Jackson, USGS, 12-9-08
0018 % Last modified: F.L. Engel, USGS 12/12/2012
0019 
0020 
0021 %% Replace bad GPS with BT
0022 
0023 for zi = 1 : z
0024     % Check for GPS data
0025     if  sum(nansum(A(zi).Nav.lat_deg)) == 0
0026         error('No GPS data detected!')
0027     end 
0028     
0029     % Prescreen GPS data (Added 4-8-10)
0030     % Eliminated uncessary "find" statement FLE 12-12-12
0031     bad_long    = A(zi).Nav.long_deg < -180 | A(zi).Nav.long_deg > 180;
0032     bad_lat     = A(zi).Nav.lat_deg < -90 | A(zi).Nav.lat_deg > 90;
0033     A(zi).Nav.long_deg(bad_long)    = NaN;
0034     A(zi).Nav.long_deg(bad_lat)     = NaN;
0035     A(zi).Nav.lat_deg(bad_long)     = NaN;
0036     A(zi).Nav.lat_deg(bad_lat)      = NaN;
0037 
0038     % Convert lat long into xUTM and yUTM
0039     [A(zi).Comp.xUTMraw,A(zi).Comp.yUTMraw,A(zi).Comp.utmzone] = ...
0040         deg2utm(A(zi).Nav.lat_deg,A(zi).Nav.long_deg);
0041     
0042     %Check if data spans two UTM zones (Added 1-16-09 P.R. Jackson)
0043     %(requires mapping toolbox)
0044     multiple_utm_zones = strmatch(A(zi).Comp.utmzone(:,1), A(zi).Comp.utmzone);
0045 %     if sum(indx) < length(A(zi).Comp.utmzone)
0046 %         disp('Caution:  Data Spans Multiple UTM Zones')  %This is not
0047 %         %working yet (1-16-09)
0048 %     end
0049     
0050 %     [latlim,lonlim] = utmzone(A(zi).Comp.utmzone(1));
0051 %     latindx = find(A(zi).Nav.lat_deg < latlim(1) | A(zi).Nav.lat_deg > latlim(2));
0052 %     lonindx = find(A(zi).Nav.lon_deg < lonlim(1) | A(zi).Nav.lon_deg > lonlim(2));
0053 %     if ~isempty(latindx) | ~isempty(lonindx)
0054 %         beep
0055 %         disp('Caution:  Data Spans Multiple UTM Zones')
0056 %     end
0057 
0058     % Check for repeat values of gps position
0059     chk        = [1; diff(A(zi).Comp.xUTMraw)];
0060     repeat_gps = chk==0;
0061     
0062     % Check for dropped ensembles
0063     dropped_ensembles = isnan(A(zi).Comp.xUTMraw) | isnan(A(zi).Comp.yUTMraw);
0064     
0065     % Check for GPS fly-aways
0066     % Delta time
0067     delta_time_elapsed = [0; diff(A(zi).Sup.timeElapsed_sec)];
0068     
0069     % The fly away filter will only consider "valid" numbers in the GPS
0070     % data. First, screen out any NaNs and replace with 0.
0071     x = A(zi).Comp.xUTMraw; %x(dropped_ensembles) = 0;
0072     y = A(zi).Comp.yUTMraw; %y(dropped_ensembles) = 0;
0073     
0074     % Compute the velocity of the location from start bank to end bank
0075     dist1   = abs([0; hypot(y(1:end-1) - y(2:end),x(1:end-1) - x(2:end))]);
0076     vel1    = dist1./delta_time_elapsed;
0077     
0078     % Compute the velocity of the location from end bank to start bank
0079     dist2   = abs([hypot(y(2:end) - y(1:end-1),x(2:end) - x(1:end-1)); 0]);
0080     vel2    = dist2./delta_time_elapsed;
0081     
0082     % Compute a threshold velocity based on the data
0083     vel_threshold   = 1.5*(nanstd([dist1;dist2]) + nanmedian([dist1;dist2]))./...
0084         (std(delta_time_elapsed) + median(delta_time_elapsed));
0085     suspect1        = (dist1>vel_threshold);
0086     suspect2        = (dist2>vel_threshold);
0087     
0088     % Logical index of fly-aways
0089     %fly_aways       = suspect1 | suspect2;
0090     fly_aways       = suspect1 & suspect2;
0091     
0092     % Create a logical index of TRUE for bad values
0093     %A(zi).Comp.bv =(isnan(A(zi).Comp.xUTMraw) + (chk==0)) > 0;
0094     A(zi).Comp.gps_reject_locations = ...
0095         bad_lat             |...
0096         bad_long            |...
0097         dropped_ensembles   |...
0098         repeat_gps          |...
0099         fly_aways;
0100     
0101     A(zi).Comp.gps_fly_aways         = fly_aways;
0102     A(zi).Comp.gps_dropped_ensembles = dropped_ensembles; 
0103     A(zi).Comp.gps_repeat_locations  = repeat_gps;    
0104                 
0105     %Identify and interpolate missing bottom track data  (Added 9-8-09,
0106     %P.R. Jackson)  Required to stop bad data return when missing GPS and
0107     %bottom track.
0108     bvbt_indx                           = find(A(zi).Nav.totDistEast == -32768);
0109     A(zi).Nav.totDistEast(bvbt_indx)    = NaN;
0110     A(zi).Nav.totDistENorth(bvbt_indx)  = NaN;
0111     gvbt_indx                           = ~isnan(A(zi).Nav.totDistEast);
0112     A(zi).Nav.totDistEast(bvbt_indx)    = interp1(A(zi).Sup.ensNo(gvbt_indx),A(zi).Nav.totDistEast(gvbt_indx),A(zi).Sup.ensNo(bvbt_indx));
0113     A(zi).Nav.totDistNorth(bvbt_indx)   = interp1(A(zi).Sup.ensNo(gvbt_indx),A(zi).Nav.totDistNorth(gvbt_indx),A(zi).Sup.ensNo(bvbt_indx));
0114 
0115     % If bad values exist replace them with Bottom Track
0116     if any(A(zi).Comp.gps_reject_locations)
0117 
0118         % Bracket bad sections, identifying the index of corresponding to
0119         % the first good data point (bg) and last good data point (ed) for
0120         % each segment
0121         [I,~]   = ind2sub(size(A(zi).Comp.gps_reject_locations),find(A(zi).Comp.gps_reject_locations==1));
0122         df      = diff(I);
0123         nbrk    = sum(df>1)+1;
0124         [I2,~]  = ind2sub(size(df),find(df>1));
0125         bg(1)   = (I(1)-1);
0126         for n = 2 : nbrk
0127             bg(n) = (I(I2(n-1)+1)-1);
0128         end
0129         for n = 1 : nbrk -1
0130             ed(n) = (I(I2(n))+1);
0131         end
0132         ed(nbrk)  = I(end)+1;
0133 
0134         % Now, determine position based on bottom track starting from the
0135         % beginning and end of bad sections then average them together
0136         
0137         % Create a canvas through preallocation
0138         A(zi).Comp.xUTMf    = NaN(length(A(zi).Comp.xUTMraw),1);
0139         A(zi).Comp.xUTMb    = NaN(length(A(zi).Comp.xUTMraw),1);
0140         A(zi).Comp.yUTMf    = NaN(length(A(zi).Comp.yUTMraw),1);
0141         A(zi).Comp.yUTMb    = NaN(length(A(zi).Comp.yUTMraw),1);
0142         BG                  = NaN(length(A(zi).Comp.yUTMraw),1);
0143         ED                  = NaN(length(A(zi).Comp.yUTMraw),1);
0144 
0145         xUTM                = NaN(length(A(zi).Comp.xUTMraw),3);
0146         yUTM                = NaN(length(A(zi).Comp.yUTMraw),3);
0147 
0148         % Screen bad locations out
0149         A(zi).Comp.xUTM     = A(zi).Comp.xUTMraw;
0150         A(zi).Comp.yUTM     = A(zi).Comp.yUTMraw;
0151         A(zi).Comp.xUTM(A(zi).Comp.gps_reject_locations)   = NaN;
0152         A(zi).Comp.yUTM(A(zi).Comp.gps_reject_locations)   = NaN;
0153 
0154         for i = 1 : nbrk
0155             for j = bg(i)+1 : ed(i)-1
0156                 if bg(i) > 0
0157                     A(zi).Comp.xUTMf(j,1)               =...
0158                         A(zi).Comp.xUTMraw(bg(i),1)     -...
0159                         A(zi).Nav.totDistEast(bg(i),1)  +...
0160                         A(zi).Nav.totDistEast(j,1);
0161                     A(zi).Comp.yUTMf(j,1)               =...
0162                         A(zi).Comp.yUTMraw(bg(i),1)     -...
0163                         A(zi).Nav.totDistNorth(bg(i),1) +...
0164                         A(zi).Nav.totDistNorth(j,1);
0165                     ED(j,1)=((j-bg(i))./(ed(i)-bg(i)));
0166                     BG(j,1)=((ed(i)-j)./(ed(i)-bg(i)));
0167                 end
0168                 if ed(i) < length(A(zi).Nav.lat_deg)
0169                     A(zi).Comp.xUTMb(j,1)               =...
0170                         A(zi).Comp.xUTMraw(ed(i),1)     -...
0171                         A(zi).Nav.totDistEast(ed(i),1)  +...
0172                         A(zi).Nav.totDistEast(j,1);
0173                     A(zi).Comp.yUTMb(j,1)               =...
0174                         A(zi).Comp.yUTMraw(ed(i),1)     -...
0175                         A(zi).Nav.totDistNorth(ed(i),1) +...
0176                         A(zi).Nav.totDistNorth(j,1);
0177                 end
0178             end
0179 
0180             if  bg(i) > 0 && ed(i) < length(A(zi).Nav.lat_deg)
0181                 xUTM(:,1)=A(zi).Comp.xUTM;
0182                 xUTM(:,2)=A(zi).Comp.xUTMf.*BG;
0183                 xUTM(:,3)=A(zi).Comp.xUTMb.*ED;
0184                 yUTM(:,1)=A(zi).Comp.yUTM;
0185                 yUTM(:,2)=A(zi).Comp.yUTMf.*BG;
0186                 yUTM(:,3)=A(zi).Comp.yUTMb.*ED;
0187                 xUTM(xUTM==0)=NaN;
0188                 yUTM(yUTM==0)=NaN;
0189                 A(zi).Comp.xUTM=nansum(xUTM,2);
0190                 A(zi).Comp.yUTM=nansum(yUTM,2);
0191 
0192                 A(zi).Comp.xUTMf=NaN(length(A(zi).Comp.xUTMraw),1);
0193                 A(zi).Comp.xUTMb=NaN(length(A(zi).Comp.xUTMraw),1);
0194                 A(zi).Comp.yUTMf=NaN(length(A(zi).Comp.yUTMraw),1);
0195                 A(zi).Comp.yUTMb=NaN(length(A(zi).Comp.yUTMraw),1);
0196                 xUTM=NaN(length(A(zi).Comp.xUTMraw),3);
0197                 yUTM=NaN(length(A(zi).Comp.yUTMraw),3);
0198             else
0199                 xUTM(:,1)=A(zi).Comp.xUTM;
0200                 xUTM(:,2)=A(zi).Comp.xUTMf;
0201                 xUTM(:,3)=A(zi).Comp.xUTMb;
0202                 yUTM(:,1)=A(zi).Comp.yUTM;
0203                 yUTM(:,2)=A(zi).Comp.yUTMf;
0204                 yUTM(:,3)=A(zi).Comp.yUTMb;
0205                 xUTM(xUTM==0)=NaN;
0206                 yUTM(yUTM==0)=NaN;
0207                 A(zi).Comp.xUTM=nanmean(xUTM,2);
0208                 A(zi).Comp.yUTM=nanmean(yUTM,2);
0209 
0210                 A(zi).Comp.xUTMf=NaN(length(A(zi).Comp.xUTMraw),1);
0211                 A(zi).Comp.xUTMb=NaN(length(A(zi).Comp.xUTMraw),1);
0212                 A(zi).Comp.yUTMf=NaN(length(A(zi).Comp.yUTMraw),1);
0213                 A(zi).Comp.yUTMb=NaN(length(A(zi).Comp.yUTMraw),1);
0214                 xUTM=NaN(length(A(zi).Comp.xUTMraw),3);
0215                 yUTM=NaN(length(A(zi).Comp.yUTMraw),3);
0216             end
0217         end
0218     else
0219         A(zi).Comp.xUTM=A(zi).Comp.xUTMraw;
0220         A(zi).Comp.yUTM=A(zi).Comp.yUTMraw;
0221     end
0222     
0223     keep A z zi
0224 %     clear I I2 J J2 bg chk df ed i j nbrk sbt xUTM yUTM n zi BG ED
0225 %     clear bad_lat bad_long repeat_gps fly_aways suspect*
0226 %     clear dist* bvbt_indx gvbt_indx multiple_utm_zones vel* delta_time_elapsed
0227 end

Generated on Thu 21-Aug-2014 10:40:31 by m2html © 2005