function gps = read_gprm_FATIMA_2022(dfl,PosLims)
%{
    reads gprm file specified by dfl and returns array with the
    following columns:
    
    1   dnum      10 Hz timestamp
    2   cog         cog, deg
    3   sog         sog, m/s
    4   Ngps        filtered N component of ship speed, m/s
    5   Egps        filtered E component of ship speed, m/s
    6   lat
    7   lon
    
    input parameters: dfl = string path to scs file
                      PosLims: [lon_lo,lon_hi,lat_lo,lat_hi] limits
    
    gps files are reformatted to gprm using a python script.
    GPRMC lines are stripped from the raw file and saved to a new gprm file.
    
    Format for $GPRMC lines in the raw file:
    0000078 $GPRMC,000000.00,A,5822.95924260,N,04515.92776422,W,1.58,223.96,151013,22.0,W,A*04
    1 2 3                      4             5 6              7 8    9      10     11
    PSD timecodes              lat             lon              spd  hdg
    
    PYTHON code:
    
    """HiWinGS PSD gps file parse: parseGpsFiles.py
    Parses data from Hemisphere gps files saved on the PSD data system.
    saves GPRMC lines with the following format to an output file.
    0000078 $GPRMC,000000.00,A,5822.95924260,N,04515.92776422,W,1.58,223.96,151013,22.0,W,A*04
    
    source file names follow the convention: gps0yydddhh_raw.txt
    output file names follow the convention: gprmyydddhh_raw.txt
    
    Using Anaconda iPython installation.
    Versions:anaconda 1.6.0, python 2.7.5
    
    Byron Blomquist, NOAA/ESRL/PSD3, Aug 2014
    """
    import os
    
    # assume current directory contains raw gps files,
    # dump list of data file names into a list variable
    files = [f for f in os.listdir('.') if f.startswith('gps')]
    
    # parse data files and save into txt files
    for file in files:
        # open output file, new name begins with gprmc
        outp = 'gprm'+file[4:]
        fgprmc = open(outp, 'w')
        print file
        fin = open(file)
        for line in fin:
            cleanLine = line.strip()
            lineOut = cleanLine + '\n'
            fields = cleanLine.split(',')
            if '$GPRMC' in fields[0]:
                fgprmc.write(lineOut)
    
        print outp+' finished'
    
    fgprmc.close()
%}

[~,fname,~] = fileparts(dfl);
YY = fname(5:6); yyyy = 2000 + str2double(YY);
DDD = fname(7:9); ddd = str2double(DDD);
HH = fname(10:11); hh = str2double(HH);

%% 10 Hz reference timestamp and averaging bins
start = datenum(yyyy-1,12,31) + ddd + hh/24;
delta = double(1.0/864000);
last = start + 35999*delta;
dnum = start:delta:last;  dnum = dnum(:); % exact 10Hz timestamp
dnum_bin = start:delta:last + double(1.0/864000);
dnum_bin = dnum_bin(:);

%% read data if file exists, if not return NaNs
if exist(dfl,'file')==2
    disp(['Reading gprm file for ',DDD,' ',HH]);
    %% read file
    flist = fopen(dfl);
    temp = {};    % empty arrays for textscan
    stx = [];

    while feof(flist)==0    % read entire file into cell array
        try
            % temp will be 1x11 cell array, cells will be ~36000x1 arrays
            temp = textscan(flist,'%2f%2f%3f %*6c %*9c %*c %f %c %f %c %f %f %f %f %*c %*4c %*[^\n]','delimiter',', ','headerlines',0,'emptyvalue',NaN);
            temp{5} = double(temp{5}); % will be ascii code for N or S
            temp{7} = double(temp{7}); % will be ascii code for E or W
            stx = [stx cell2mat(temp)'];  % 11x~36000 array
        catch
            for ii=1:11      % length of last cell reflects missing values from any/all fields
                 if length(temp{1,ii})~=length(temp{1,11})
                     temp{1,ii}(length(temp{1,11})+1)=[];  %truncate length of all cells to length(temp{1,11})
                 end
            end
            if ~isempty(temp{1,1})  % save good values from last read, if any
                temp{5} = double(temp{5});
                temp{7} = double(temp{7});
                stx = [stx cell2mat(temp)'];
            end
        end
    end
    fclose(flist);
    stx = stx';
    dmun0 = start + stx(:,1)/1440 +(stx(:,2)+stx(:,3)/1000)/86400;

    %% format variables
    glt = stx(:,4);
    gpslat_sign = stx(:,5);             % ascii code for "N" or "S"
    gln = stx(:,6);
    gpslon_sign = stx(:,7);             % ascii code for "W" or "E"
    gpsspeed = stx(:,8);                % in kts
    gpsspeed(gpsspeed>20) = NaN;        % kill wierd values
    gpsspeed(gpsspeed<0) = NaN; 
    gpshead = stx(:,9);                 % in deg
    gpshead(gpshead>360) = NaN;         % kill wierd values
    gpshead(gpshead<0) = NaN;
    gpshead = unwrap(gpshead*pi/180);   % unwrapped in radians for now...

    bb = floor(glt/100);                % convert lat format to decimal degrees
    glt = (glt-bb*100)/60;
    glt = glt+bb;
    ii = find(gpslat_sign==83); glt(ii) = -glt(ii);

    bbn = floor(gln/100);               % convert lon format to decimal degrees
    gln = (gln-bbn*100)/60;
    gln = gln+bbn;
    ii = find(gpslon_sign==87);         % W longitudes
    gln(ii) = -gln(ii) + 360;           % reformat to E longitude, 0-360 deg

    % discard unreasonable lat/lon values
    gln(gln<PosLims(1)) = NaN;
    gln(gln>PosLims(2)) = NaN;
    glt(glt<PosLims(3)) = NaN;
    glt(glt>PosLims(4)) = NaN;

    %% despike variables
    glt = despike2(glt);
    gln = despike2(gln);
    gpshead = despike2(gpshead);
    gpsspeed = despike2(gpsspeed);

    % interp to dnum grid
    [~,zz,~] = unique(dmun0);
    gpslon1 = interp1(dmun0(zz),gln(zz),dnum,'linear');
    gpslat1 = interp1(dmun0(zz),glt(zz),dnum,'linear');
    gpshead1 = interp1(dmun0(zz),gpshead(zz),dnum,'linear');
    gpsspeed1 = interp1(dmun0(zz),gpsspeed(zz),dnum,'linear');

    % look for spikes again
    ii = abs(diff(gpslon1))>0.01; gpslon1(ii) = NaN;
    ii = abs(diff(gpslat1))>0.01; gpslat1(ii) = NaN;

    % replace NaNs with nearest neighbor
    % doesn't 
    gpslon1 = replace_NaN_nearest_neighbor(gpslon1);
    gpslat1 = replace_NaN_nearest_neighbor(gpslat1);
    gpsspeed1 = replace_NaN_nearest_neighbor(gpsspeed1);
    gpshead1 = replace_NaN_nearest_neighbor(gpshead1);

    cog = gpshead1.*180/pi;    % back to degrees
    cog = mod(cog+360,360);    % be sure range is 0-360
    sog = gpsspeed1*0.5144;    % in m/s

    %% compute ship N & E velocity components in m/s
    Ngps = sog.*cos(cog./180.*pi); % 1x36000 arrays
    Egps = sog.*sin(cog./180.*pi);

    % filter ship velocities unless excessive NaNs are present
    if all(isfinite(Ngps))
        fs = 10;         % GPS output @ 10 Hz
        fc = 1/200;      % Cutoff frequency
        wp = 1/25/(fs/2);
        ws = 0.8*wp;
        % [n,wn]=buttord(wp,ws,10,25);
        [n,wn] = buttord(wp,ws,3,7);
        [bhi,ahi] = butter(n,wn,'low');
        Ngps = filtfilt(bhi,ahi,Ngps);
        Egps = filtfilt(bhi,ahi,Egps);
    end

    %% output
    gps = [dnum cog sog Ngps Egps gpslat1 gpslon1];
else
    gps = NaN(36000,7);
    gps(:,1) = dnum;
end

% reformat to table for output
vars = {'dnum','cog','sog','Ngps','Egps','lat','lon'};
units = {'none','deg','m/s','m/s','m/s','degN','degE'};
gps = array2table(gps,'VariableNames',vars);
gps.Properties.VariableUnits = units;

end
