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 (filepath, tabular, convert_mjd)
function varargout = read_teqc_compact_each (filepath, tabular, convert_mjd)
    if (nargin < 2) || isempty(tabular),  tabular = true;  end
    if (nargin < 3),  convert_mjd = [];  end
    myrepmat = @(x, siz) cast(x .* ones(siz), class(x));  % faster.
    [convert_mjd, epoch_digits] = check_mydate (convert_mjd);
    
    %% 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 = cellfun(@(x)x(1:end-1), data2, 'UniformOutput',false);
    
    %%
    assert(strstarti('COMPACT', data2{1}))
    version = 1;
    if strstarti('COMPACT2', data2{1}),  version = 2;  end

    %%
    if (version == 1)
        header = data2(1:4);
        data2(1:4) = [];
        assert(strstarti('SVS', header{2}))
        header(2) = [];  % make it similar to version 2.
    else
        % (line "SVS" is absent.)
        header = data2(1:3);
        data2(1:3) = [];
    end
      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_mjd
        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);
    obs = cellfun(@(x) sscanf(x, '%11f%*c')', obs_str, 'UniformOutput',false);
    if (version == 1)
        sat     = cellfun(@(x) sscanf(x, '%d')', sat_str, 'UniformOutput',false);
        num_sat = cellfun(@(x)x(1), sat);
        prn     = cellfun(@(x)x(2:end), sat, 'UniformOutput',false);
        sys     = cellfun(@(x)myrepmat('G', size(x)), prn, 'UniformOutput',false);
    else
        [num_sat_cell, ignore, ignore, ind_cell] = cellfun(@(x) sscanf(x, '%d', 1), sat_str, ...
            'UniformOutput',false); %#ok<ASGLU>
        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 = cellfun(@(sat_str_i,ind_i) sat_str_i(ind_i:end), ...
            sat_str, ind_cell, 'UniformOutput',false);
        temp = cellfun(@(sat_str_i,num_sat_i) sscanf(sat_str_i, ' %c%u', max(0,num_sat_i*2))', ...
            sat_str, num_sat_cell, 'UniformOutput',false);
        prn = cellfun(@(x) x(2:2:end), temp, 'UniformOutput',false);
        sys = cellfun(@(x) x(1:2:end), temp, 'UniformOutput',false);
        sys = cellfun(@char, sys, 'UniformOutput',false);
    end
    num_sat2 = cellfun(@(x)size(x,2), prn);
    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)))

    %%
    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);
    clear num_sat_cell

    %%
    time = dt.*(0:numel(obs)-1)';
    if convert_mjd
        epoch = epoch0 + mydateseci(time);
        % otherwise epoch intersection will fail:
        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 = cellfun(@(epoch_i,num_sat_i) myrepmat(epoch_i, [num_sat_i 1]), ...
            num2cell(epoch), num2cell(num_sat), 'UniformOutput',false);        
        epoch = cell2mat(temp);
    end
    
    %%    
    [varargout{1:nargout}] = read_teqc_compact_out (obs, epoch, prn, sys);

    %%
    % File format reference:
    % <http://ls.unavco.org/pipermail/teqc/2009/000827.html>
    % <http://ls.unavco.org/pipermail/teqc/2007/000566.html>
end

function [convert_mjd, epoch_digits] = check_mydate (convert_mjd)
    epoch_digits = NaN;
    convert_mjd_input = convert_mjd;
    if isempty(convert_mjd),  convert_mjd = true;  end
    if ~convert_mjd,  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_mjd_input)  % && (convert_mjd_input == true)
            error('MATLAB:read_teqc_compact:MJDnotPossible', ...
                'Date functions unavailable.');
        end
        %warning('MATLAB:read_teqc_compact:MJD', 'Output epoch is in MJD.');
        convert_mjd = false;
    end
end

Contact us