# 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:

- VMT_PreProcess This function is a driver to preprocess the transects data. Several

## 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