image thumbnail

read_teqc_compact

by

 

17 Apr 2011 (Updated )

Reads teqc's QC compact files more rapidly -- without for loops.

read_teqc_compact_each
function varargout = read_teqc_compact_each (filepath, tabular, convert_epoch)
    %disp(filepath)  % DEBUG
    if (nargin < 2) || isempty(tabular),  tabular = true;  end
    if (nargin < 3),  convert_epoch = [];  end
    myrepmat = @(x, siz) cast(x .* ones(siz), class(x));  % faster?
    myrepmat = @repmat;  % DEBUG
    [convert_epoch, epoch_digits] = check_mydate (convert_epoch);
    cellfun2 = @(varargin) cellfun(varargin{:}, 'UniformOutput',false);
    rowvec = @(x) x(:).';
    
    %% read the whole file content at once:
    [fid, temp_msg] = fopen (filepath, 'rt');
      %[filename, permission, machineformat, encoding] = fopen(fid)  % DEBUG
    if (fid == -1)
        error('MATLAB:read_teqc_compact:badFile', ...
            'Could not open file "%s":\n%s', filepath, temp_msg);
    end
    clear temp_msg
    data = fread(fid, '*char')';
    fclose (fid);
    if strcmp(data(end), sprintf('\n')),  data(end) = [];  end

    %% "if the epoch lines lists 0 SVs, then there is no following data line"
    % hint: try inserting fake empty data line, before doing the parsing:
    %data = strrep(data, sprintf('\n 0\n'), sprintf('\n 0\n\n'));  % WRONG!
    data = strrep(data, sprintf('\n 0'), sprintf('\n 0\n'));

    %% split up lines:
    data = reshape(data, 1, []);
    temp = find(data == sprintf('\n'));
    temp2 = diff([0; temp(:); numel(data)]);
    data2 = mat2cell(data, 1, temp2).';
    % discard newline character:
    % (unnecessary)
    %data2 = cellfun2(@(x)x(1:end-1), data2);
    
    %%
    switch upper(strtrim(data2{1}))
    case 'COMPACT',   version = 1;
    case 'COMPACT2',  version = 2;
    case 'COMPACT3',  version = 3;
    otherwise, error('MATLAB:read_teqc_compact:badFile', ...
        'Bad first header line in file "%s":\n%s', filepath);
    end
            
    %%
    switch version
    case 1
        header = data2(1:4);
        data2(1:4) = [];
        assert(strstarti('SVS', header{2}));
        header(2) = [];  % make it similar to version 2.
    case 2
        % (line "SVS" is absent.)
        header = data2(1:3);
        data2(1:3) = [];
    case 3
        % (line "SVS" and line "T_SAMP" are absent.)
        header = data2(1:2);
        data2(1:2) = [];
    end
    
    if (version == 1 || version == 2)
        assert(iseven(size(data2,1)));
        %celldisp(header)  % DEBUG
  
        %%
        assert(strstarti('T_SAMP', header{2}));
        dt = sscanf(header{2}, 'T_SAMP %f', [1 1]);
        
        %%
        assert(strstarti('START_TIME_', header{3}))
        epoch0 = sscanf(header{3}, 'START_TIME_MJ%*c %f', [1 1]);
        % "the start time in modified Julian days is now labeled 'START_TIME_MJD'"
        % (in version 2);  in version 1 it was 'START_TIME_MJL'.
        if convert_epoch
            epoch0 = mydatemjdi(epoch0);
            %epoch0 = roundn(epoch0, epoch_digits);  (not enough; needs to do for all epochs)
        end

        %%
        sat_str = data2(1:2:end);
        obs_str = data2(2:2:end);
    else % (version == 3)
        assert(strstarti('GPS_START_TIME', header{2}));
        epoch0 = sscanf(header{2}, 'GPS_START_TIME %d %d %d %d %f %f', [1 6]);
        epoch0 = double(epoch0);
        if convert_epoch
            epoch0 = mydatenum(epoch0);
        else
            epoch0 = 0;
        end
        rec_str = data2(1:2:end);
        obs_str = data2(2:2:end);
        time_str = cellfun2(@(x) x(1:11), rec_str);
        sat_str = cellfun2(@(x) x(13:end), rec_str);
    end
    
    obs = cellfun2(@(x) sscanf(x, '%11f%*c')', obs_str);
    if (version == 1)
        sat     = cellfun2(@(x) sscanf(x, '%d')', sat_str);
        %num_sat = cellfun2(@(x)x(1), sat);  % WRONG!
        num_sat = cellfun(@(x)x(1), sat);
        prn     = cellfun2(@(x)x(2:end), sat);
        num_sat2 = cellfun(@(x) size(x,2), prn);
        sys     = cellfun2(@transpose, mat2cell(myrepmat('G', [sum(num_sat2) 1]), num_sat2));
        %sys     = cellfun2(@(x) myrepmat('G', size(x)), prn);  % slower
    else % (version == 2 || version == 3)
        [num_sat_cell, ~, ~, ind_cell] = cellfun2(@(x) sscanf(x, '%d', 1), sat_str);
        num_sat = cell2mat(num_sat_cell);
        %ind = cell2mat(ind_cell);
        %fnc = @(i) sscanf(sat_str{i}(ind(i):end), ' %c%u', max(0,num_sat(i)*2))';
        %sat_str = arrayfun(fnc, (1:numel(sat_str))', 'UniformOutput',false);
        sat_str = cellfun2(@(sat_str_i,ind_i) sat_str_i(ind_i:end), sat_str, ind_cell);
        temp = cellfun2(@(sat_str_i,num_sat_i) sscanf(sat_str_i, ' %c%u', max(0,num_sat_i*2))', ...
            sat_str, num_sat_cell);
        prn = cellfun2(@(x) x(2:2:end), temp);
        sys = cellfun2(@(x) x(1:2:end), temp);
        sys = cellfun2(@char, sys);
        num_sat2 = cellfun(@(x)size(x,2), prn);
        clear num_sat_cell
    end
    idx = (num_sat <= 0);  % = ( (num_sat == -1) | (num_sat == 0) );
    assert(all(num_sat2(idx)==0));
    %[num_sat2(~idx), num_sat(~idx)]  % DEBUG
    assert(isequal(num_sat2(~idx), num_sat(~idx)));

    %%
    num_sat3 = cellfun(@numel, obs);
    idx3 = (num_sat ~= num_sat3) & ~idx;
    if any(idx3)
        warning('MATLAB:read_teqc_compact:unSync', ...
            ['Spurious observations detected; discarding:'...
             '\n(epochs: %s)\n\t%s'], str2list(time_str(idx3)), filepath);
        for ind=rowvec(find(idx3))
            obs{ind}(num_sat(ind)+1:end) = [];
        end
    end
    
    %%
    idx = (num_sat == 0);
      assert(all(cellfun(@isempty, obs(idx))));
    obs(idx) = {[]};  % it was '' now is []; otherwise cell2mat yields error 'MATLAB:cell2mat:MixedDataTypes'.
    prn(idx) = {[]};
      
    %% "a -1 for the number of SVs means the SV list is the same as the previous epoch":
    temp = find(num_sat~=-1);
    temp2 = diff([temp; numel(num_sat)+1]);
    temp3 = cell2mat(arrayfun(@(x,y) myrepmat(x, [y,1]), temp, temp2, 'UniformOutput',false));
      assert(numel(temp3) == numel(num_sat))
    prn = prn(temp3);
    sys = sys(temp3);
    num_sat = num_sat(temp3);

    %%
    if (version == 1 || version == 2)
        time = dt.*(0:numel(obs)-1)';
    else % (version == 3)
        time = str2double(time_str);
    end
    if convert_epoch
        epoch = epoch0 + mydateseci(time);
        % necessary to prevent failed epoch intersection:
        epoch = roundn(epoch, epoch_digits);
    else
        epoch = epoch0 + time ./ (24*60^2);
    end

    %%
    if tabular
        obs = cell2mat(obs')';
        prn = cell2mat(prn')';
        sys = cell2mat(sys')';
        %%%temp = arrayfun(@(epoch_i,num_sat_i) myrepmat(epoch_i, [num_sat_i 1]), ...
        %%%    epoch, num_sat, 'UniformOutput',false);
        %%temp = cellfun2(@(epoch_i,num_sat_i) myrepmat(epoch_i, [num_sat_i 1]), ...
        %%    num2cell(epoch), num2cell(num_sat));
        %num_epochs = numel(epoch);
        %temp = cell(num_epochs,1);
        %for i=1:num_epochs,  temp{i} = myrepmat(epoch(i), [num_sat(i) 1]);  end
        %epoch = cell2mat(temp);
        num_epochs = numel(epoch);
        max_num_sat = max(num_sat);
        temp = myrepmat(rowvec(epoch), [max_num_sat 1]);
        for i=1:num_epochs,  temp(num_sat(i)+1:max_num_sat,i) = NaN;  end
        epoch = temp(:);
        epoch(isnan(epoch)) = [];
    end
       
    %%    
    [varargout{1:nargout}] = read_teqc_compact_out (obs, epoch, prn, sys);
    
    %%
    % File format reference:
    % <http://postal.unavco.org/pipermail/teqc/2013/001594.html>
    % <http://ls.unavco.org/pipermail/teqc/2009/000827.html>
    % <http://ls.unavco.org/pipermail/teqc/2007/000566.html>
end

function [convert_epoch, epoch_digits] = check_mydate (convert_epoch)
    epoch_digits = NaN;
    convert_epoch_input = convert_epoch;
    if isempty(convert_epoch),  convert_epoch = true;  end
    if ~convert_epoch,  return;  end
    try
        epoch_precision = mydatemjdi(0.000000001/2) - mydatemjdi(0);
        epoch_digits = log10(epoch_precision);
        epoch_digits = fix(epoch_digits);
    catch ME
        if ~strcmp(ME.identifier, 'MATLAB:UndefinedFunction'),  rethrow(ME);  end
        if ~isempty(convert_epoch_input)  % && (convert_epoch_input == true)
            error('MATLAB:read_teqc_compact:MJDnotPossible', ...
                'Date functions unavailable.');
        end
        %warning('MATLAB:read_teqc_compact:MJD', 'Output epoch is in MJD.');
        convert_epoch = false;
    end
end

Contact us