Fun_LoadWaveStation.m

% =========================================================================
% Load station information and three-component waveform files for
% finite-fault inversion.
%
% Input:
%   outStaFolder  : Folder containing StationInfo.dat
%   outWaveFolder : Folder containing station waveform files (*.dat)
%   targetRate    : Desired sampling rate [Hz]
%
% Output:
%   ob            : Observed acceleration data
%                   Size = [numberOfSamples, 3, numberOfStations]
%                   Component order = [EW, NS, UD]
%
%   loca          : Station coordinates [latitude, longitude]
%                   Size = [numberOfStations, 2]
%
%   tstart        : Station record beginning relative to origin time [s]
%
%   StaName       : Station names as a char array
%
%   rate0         : Original sampling rate in StaInfo.txt [Hz]
%
% Important:
%   This function reads:
%       OutSta/StaInfo.txt
%       OutWave/<StationName>.dat
%
%   Each waveform file must contain three columns:
%       EW  NS  UD
%
%   The waveform unit is assumed to be m/s^2.
% =========================================================================
function [ob, loca, tstart, StaName, rate0] = ...
    Fun_LoadWaveStation(outStaFolder, outWaveFolder, targetRate)
stationInfoFile = fullfile(outStaFolder, 'StaInfo.txt');
if ~exist(stationInfoFile, 'file')
    error('StaInfo.txt was not found: %s', stationInfoFile);
end
if ~isfolder(outWaveFolder)
    error('Waveform folder does not exist: %s', outWaveFolder);
end

%% ------------------------------------------------------------------------
% Read StaInfo.txt
% -------------------------------------------------------------------------
fid = fopen(stationInfoFile, 'r');
if fid < 0
    error('Unable to open StaInfo.txt.');
end
numberOfStationsInHeader = [];
dataType = [];
samplingInterval = [];

foundStationHeader = false;

while ~feof(fid)

    line = fgetl(fid);

    if ~ischar(line)
        break;
    end

    % Read the sampling interval line after its header line.
    if contains(line, '# Nsta Wave-Type(0/1/2=dis/vel/acc) Time-interval [s]')

        parameterLine = fgetl(fid);
        values = sscanf(parameterLine, '%f');

        if numel(values) >= 3
            numberOfStationsInHeader = values(1);
            dataType = values(2);
            samplingInterval = values(3);
        end
    end

    % Locate the station table header.
    if contains(line, '# station latitude longitude Tstart Nsample')
        foundStationHeader = true;
        break;
    end
end

if ~foundStationHeader
    fclose(fid);
    error(['Cannot find the station-table header in StaInfo.txt. ', ...
           'Expected: # station latitude longitude Tstart Nsample']);
end
StaInfo = textscan(fid, '%s %f %f %f %f','MultipleDelimsAsOne', true);
fclose(fid);

StaName = char(StaInfo{1});
stationLatitude = StaInfo{2};
stationLongitude = StaInfo{3};
tstart = StaInfo{4};
numberOfSamplesInInfo = StaInfo{5};
loca = [stationLatitude, stationLongitude];

if isempty(samplingInterval) || samplingInterval <= 0
    error('Invalid sampling interval in StaInfo.txt.');
end
rate0 = 1 / samplingInterval;
numberOfStations = size(loca, 1);
if ~isempty(numberOfStationsInHeader) && ...
        numberOfStations ~= numberOfStationsInHeader
    warning(['Station number mismatch: zStationInfo header says %d, ', ...
             'but the station table contains %d stations.'], ...
             numberOfStationsInHeader, numberOfStations);
end
if ~isempty(dataType) && dataType ~= 2
    warning(['The data type in StaInfo.txt is %d. ', ...
             'This script assumes type 2 means acceleration.'], dataType);
end


%% ------------------------------------------------------------------------
% Find waveform files in OutWave
% -------------------------------------------------------------------------
waveFiles = dir(fullfile(outWaveFolder, '*.dat'));
if isempty(waveFiles)
    error('No waveform files (*.dat) were found in: %s', outWaveFolder);
end
waveFileNames = string({waveFiles.name})';
waveStationNames = regexprep(waveFileNames, '\.dat$', '');
stationNameString = string(strtrim(cellstr(StaName)));
% Match waveform files according to station name, not directory order.
[isFound, waveIndex] = ismember( ...
    lower(stationNameString), lower(waveStationNames));
if any(~isFound)

    missingNames = stationNameString(~isFound);

    error(['The following stations exist in StaInfo.txt but their ', ...
           'waveform files are missing from OutWave:\n%s'], ...
           strjoin(missingNames, newline));
end

waveFiles = waveFiles(waveIndex);

if numel(waveFiles) ~= numberOfStations
    error('The number of waveform files does not match the station table.');
end

%% ------------------------------------------------------------------------
% Read, resample, and align waveforms
% -------------------------------------------------------------------------
% Calculate rational resampling factors.
[p, q] = rat(targetRate / rate0, 1e-12);
traceCell = cell(numberOfStations, 1);
sampleOffset = round(tstart * targetRate);
for iStation = 1:numberOfStations
    waveFile = fullfile(outWaveFolder, waveFiles(iStation).name);
    waveform = load(waveFile);
    if size(waveform, 2) < 3
        error(['Waveform file %s has fewer than three columns. ', ...
               'Expected EW / NS / UD.'], waveFile);
    end
    % Retain EW, NS, and UD only.
    waveform = waveform(:, 1:3);

    if any(~isfinite(waveform(:)))
        error('NaN or Inf values were found in waveform file: %s', waveFile);
    end
    % Resample all three components.
    waveform = resample(waveform, p, q);

    % Align waveform beginning to earthquake origin time.
    if sampleOffset(iStation) >= 0
        waveform = [zeros(sampleOffset(iStation), 3); waveform];
    else
        firstSample = 1 - sampleOffset(iStation);
        if firstSample > size(waveform, 1)
            warning(['Station %s begins after all available samples ', ...
                     'would be removed. The station is set to zero.'], ...
                     strtrim(StaName(iStation, :)));

            waveform = zeros(1, 3);

        else
            waveform = waveform(firstSample:end, :);
        end
    end

    traceCell{iStation} = waveform;

    % This check is only informative; it does not force truncation.
    expectedSamples = round(numberOfSamplesInInfo(iStation) * targetRate / rate0);

    if abs(size(waveform, 1) - expectedSamples - sampleOffset(iStation)) > 2
        warning(['Station %s has an unexpected waveform length after ', ...
                 'resampling/alignment.'], ...
                 strtrim(StaName(iStation, :)));
    end
end

%% ------------------------------------------------------------------------
% Build ob = [time, component, station]
%
% The final length follows the original script logic:
% use the longest available waveform and zero-pad shorter records.
% -------------------------------------------------------------------------
traceLengths = cellfun(@(x) size(x, 1), traceCell);
maximumSamples = max(traceLengths);
ob = zeros(maximumSamples, 3, numberOfStations);
for iStation = 1:numberOfStations
    waveform = traceCell{iStation};
    ob(1:size(waveform, 1), :, iStation) = waveform;
end


fprintf('[Loading] Waveforms loaded successfully\n');
fprintf('   Observed data: %dx%dx%d\n', size(ob, 1), size(ob, 2), size(ob, 3));
fprintf('   Origin Sampling: %.1f Hz\n', rate0);
fprintf('   Target Sampling: %.1f Hz\n', targetRate);
fprintf('[Loading] Stations: %d\n\n', numberOfStations);

end