Function working with helping function

Hello

I got quite a long code and I might be doing it in a slightly complicated way but I'm not sure what to improve, either way it doesn't work at the moment and I'm getting an error.

Generally - my main goal is to find drops of 2400 in value from the result_matrix, then linking these drops with the assumption that they happen next to each other (geographically i.e., one grid cell next to another) over a timespan of consecutive X days - in my code it's 7, but it's irrelevant. Then store them in another matrix. Eventually I will want to overlay the resulting 'tracks' on the map, but I'll try to figure it out myself. I'm attaching an example here (in my case it'd be much less complex).

That is my code:

% Constants
THRESHOLD_PRESSURE_DROP = 2400; % Threshold for pressure drop (in hPa*100)
MAX_DISTANCE = 1; % Maximum distance for linking (in grid cells)
MIN_TRACK_LENGTH = 7; % Minimum number of time steps for a valid track
MAX_TRACK_GAP = 1; % Maximum gap in time steps for track continuity

try
    % Load the result_matrix
    if ~exist('result_matrix', 'var')
    end

    % Calculate anomalies
    mean_pressure = mean(result_matrix, 3);
    pressure_anomalies = result_matrix - mean_pressure;

    % Tracks matrix
    storm_tracks = struct('history', [], 'current', []);

    % Iterate over time 
    for t = 1:size(result_matrix, 3)
        current_centers = detect_storm_centers(pressure_anomalies(:, :, t), THRESHOLD_PRESSURE_DROP);

        % Link centers with existing tracks or start new ones
        for i = 1:length(current_centers)
            [linked, track_index] = link_center_to_tracks(current_centers(i), storm_tracks, t, MAX_DISTANCE, MAX_TRACK_GAP);
            
            if linked
                % Update current center and history
                storm_tracks(track_index).current = current_centers(i);
                storm_tracks(track_index).history(end + 1) = t;
            else
                % Start a new track with current center and history initialized with t
                new_track = struct('history', [t], 'current', current_centers(i));
                storm_tracks(end + 1) = new_track;
            end
        end
    end

    % Filter out short and dissipated tracks
    valid_tracks = storm_tracks(arrayfun(@(x) length(x.history) >= MIN_TRACK_LENGTH, storm_tracks));


catch exception
    disp(['An error occurred: ' exception.message]);
end

% Helping functions; detect_storm_centers, link_center_to_tracks

function centers = detect_storm_centers(pressure_anomalies, threshold)
    [rows, cols] = find(pressure_anomalies < -threshold);
    centers = struct('row', num2cell(rows), 'col', num2cell(cols));
end

function [linked, track_index] = link_center_to_tracks(current_center, storm_tracks, current_time, max_distance, max_track_gap)
    linked = false;
    track_index = [];
    
    for i = 1:length(storm_tracks)
        last_time = storm_tracks(i).history(end);
        
        if abs(current_time - last_time) <= max_track_gap
            distance = sqrt((current_center.row - storm_tracks(i).current.row)^2 + ...
                            (current_center.col - storm_tracks(i).current.col)^2);
            
            if distance <= max_distance
                linked = true;
                track_index = i;
                break;
            end
        end
    end
end

It seems to go wrong when it comes to a point of using link_center_to_tracks function, but I'm not sure what exactly I'm doing wrong and the error messages don't help. I'm also attaching the data file, and the code I use to load the data if anyone would be willing to help me out with that.

clear all;

close all;

clc;

% Constants

START_ROW = 5;

END_ROW = 13;

START_COL = 11;

END_COL = 21;

NUM_DATES = 365;

% Read data

data = readmatrix('data.dat.txt');

% Generate list of row numbers for the dates

date_rows = 1:17:size(data, 1);

% Initialize a 3D matrix

result_matrix = zeros(END_ROW - START_ROW + 1, END_COL - START_COL + 1, NUM_DATES);

for i = 1:length(date_rows)

% Extract data using the defined constants

data_values = data(date_rows(i) + START_ROW - 1 : date_rows(i) + END_ROW - 1, START_COL:END_COL);

% Store data in the result_matrix

result_matrix(:, :, i) = data_values;

end

% Display a slice of the result_matrix for the first three dates

for i = 1:3

disp(['Data for Date ' num2str(i)]);

disp(result_matrix(:, :, i));

disp('------------------------');

end

I appreciate it's a lot. Thanks for your time.

Answer 1

[Voss]

This line

storm_tracks = struct('history', [], 'current', [])

storm_tracks = struct with fields:
    history: []
    current: []

makes storm_tracks a scalar struct

numel(storm_tracks) % storm_tracks is scalar (1 element)

ans = 1

with each field an empty array, so that later, in link_center_to_tracks, when you loop over the elements of storm_tracks:

for i = 1:length(storm_tracks)
    last_time = storm_tracks(i).history(end);
    % ...
end

you get the error because storm_tracks(1).history is empty. That is, there is no 'end' element in storm_tracks(1).history.

Instead, initialize storm_tracks to be an empty struct array - something like this:

storm_tracks = struct('history', {}, 'current', {})

storm_tracks = 

  0×0 empty struct array with fields:

    history
    current

numel(storm_tracks) % storm_tracks is empty (0 elements)

ans = 0

That way, when you loop over storm_tracks, the loop executes zero times because storm_tracks has zero elements, and no error occurs. linked will be false and will be returned to the main script, so a new element will be appended to storm_tracks, which has non-empty history and current fields, as intended, and the code will run the rest of the way.

Here it is with that change to how storm_tracks is initialized. You can see you end up with 37 valid_tracks.

clear all;
close all;
clc;

% Constants
START_ROW = 5;
END_ROW = 13;
START_COL = 11;
END_COL = 21;
NUM_DATES = 365;

% Read data
data = readmatrix('data.dat.txt');

% Generate list of row numbers for the dates
date_rows = 1:17:size(data, 1);

% Initialize a 3D matrix
result_matrix = zeros(END_ROW - START_ROW + 1, END_COL - START_COL + 1, NUM_DATES);

for i = 1:length(date_rows)
    % Extract data using the defined constants
    data_values = data(date_rows(i) + START_ROW - 1 : date_rows(i) + END_ROW - 1, START_COL:END_COL);
    
    % Store data in the result_matrix
    result_matrix(:, :, i) = data_values;
end

% Display a slice of the result_matrix for the first three dates
% for i = 1:3
    % disp(['Data for Date ' num2str(i)]);
    % disp(result_matrix(:, :, i));
    % disp('------------------------');
% end

% Constants
THRESHOLD_PRESSURE_DROP = 2400; % Threshold for pressure drop (in hPa*100)
MAX_DISTANCE = 1; % Maximum distance for linking (in grid cells)
MIN_TRACK_LENGTH = 7; % Minimum number of time steps for a valid track
MAX_TRACK_GAP = 1; % Maximum gap in time steps for track continuity

try
    % Load the result_matrix
    if ~exist('result_matrix', 'var')
    end
    
    % Calculate anomalies
    mean_pressure = mean(result_matrix, 3);
    pressure_anomalies = result_matrix - mean_pressure;

    % Tracks matrix
    % storm_tracks = struct('history', [], 'current', []);
    storm_tracks = struct('history', {}, 'current', {});

    % Iterate over time 
    for t = 1:size(result_matrix, 3)
        current_centers = detect_storm_centers(pressure_anomalies(:, :, t), THRESHOLD_PRESSURE_DROP);
        
        % Link centers with existing tracks or start new ones
        for i = 1:length(current_centers)
            [linked, track_index] = link_center_to_tracks(current_centers(i), storm_tracks, t, MAX_DISTANCE, MAX_TRACK_GAP);
            
            if linked
                % Update current center and history
                storm_tracks(track_index).current = current_centers(i);
                storm_tracks(track_index).history(end + 1) = t;
            else
                % Start a new track with current center and history initialized with t
                new_track = struct('history', [t], 'current', current_centers(i));
                storm_tracks(end + 1) = new_track;
            end
        end
    end

    % Filter out short and dissipated tracks
    valid_tracks = storm_tracks(arrayfun(@(x) length(x.history) >= MIN_TRACK_LENGTH, storm_tracks))


catch exception
    disp(['An error occurred: ' exception.message]);
end

valid_tracks = 1×37 struct array with fields:
    history
    current

% Helping functions; detect_storm_centers, link_center_to_tracks

function centers = detect_storm_centers(pressure_anomalies, threshold)
    [rows, cols] = find(pressure_anomalies < -threshold);
    centers = struct('row', num2cell(rows), 'col', num2cell(cols));
end

function [linked, track_index] = link_center_to_tracks(current_center, storm_tracks, current_time, max_distance, max_track_gap)
    linked = false;
    track_index = [];
    
    for i = 1:length(storm_tracks)
        last_time = storm_tracks(i).history(end);
        
        if abs(current_time - last_time) <= max_track_gap
            distance = sqrt((current_center.row - storm_tracks(i).current.row)^2 + ...
                            (current_center.col - storm_tracks(i).current.col)^2);
            
            if distance <= max_distance
                linked = true;
                track_index = i;
                break;
            end
        end
    end
end

Note: I did not examine the logic of the code beyond the cause of and resolution to the error, so I make no claim that the code does what it's supposed to do, but now it should run and you can determine whether the results are correct.

[OcDrive]

Perfect, thank you Voss.

[Voss]

You're welcome OcDrive.