Matlabify - River Module

mhkit/river/IO/read_usgs_file.m

@@ -1,56 +1,112 @@
-function datast=read_usgs_file(file_name)
+function dataset=read_usgs_file(file_name)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Reads a USGS JSON data file (from https://waterdata.usgs.gov/nwis)
-%     into a structure
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Reads a USGS JSON data file (from https://waterdata.usgs.gov/nwis)
+% into a structure
 %
 % Parameters
-% ----------
-%     file_name : str
-%         Name of USGS JSON data file
+% ------------
+% file_name : string
+%   Name of USGS JSON data file
 %
 % Returns
-% -------
-%     datast : structure
-%
-%
-%         datast.Data: named according to the parameter's variable description
+% ---------
+% dataset : struct
+%   USGS data, one field per parameter
+%     dataset.(parameter) : column vector
+%       Data named according to the parameter's variable description
+%       (e.g. dataset.Discharge)
+%     dataset.time : row vector [s]
+%       Epoch time
+%     dataset.units : struct
+%       Units for each parameter (e.g. dataset.units.Discharge)
 %
-%         datast.time: epoch time [s]
-%
-%         datast.units: units for each parameter
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
+arguments (Input)
+    file_name {mustBeTextScalar, mustBeFile}
+end
 
-datapd=py.mhkit.river.io.usgs.read_usgs_file(file_name);
+arguments (Output)
+    dataset struct
+end
 
+text = jsondecode(fileread(file_name));
 
-xx=cell(datapd.axes);
-v=xx{2};
+time_series = text.value.timeSeries;
+num_series = numel(time_series);
 
+names = strings(num_series, 1);
+units = strings(num_series, 1);
+series_time = cell(num_series, 1);
+series_value = cell(num_series, 1);
+keep = false(num_series, 1);
 
-vv=cell(py.list(py.numpy.nditer(v.values,pyargs("flags",{"refs_ok"}))));
+for i = 1:num_series
+    try
+        series = time_series(i);
 
-vals=double(py.array.array('d',py.numpy.nditer(datapd.values,pyargs("flags",{"refs_ok"}))));
-sha=cell(datapd.values.shape);
-x=int64(sha{1,1});
-y=int64(sha{1,2});
+        % USGS variableDescription is formatted as "<name>, <units>"
+        description = string(series.variable.variableDescription);
+        name_parts = split(description, ",");
+        names(i) = name_parts(1);
+        if numel(name_parts) > 1
+            units(i) = name_parts(2);
+        end
 
-vals=reshape(vals,[x,y]);
-si=size(vals);
+        % Only the first entry of "values" is used, matching the Python
+        % implementation's use of values[0]["value"]
+        values = series.values(1).value;
 
-for i=1:si(2)
-    test=string(py.str(vv{i}));
-    newname=split(test,",");
-    datast.(newname(1))=vals(:,i);
-    datast.units.(newname(1))=newname(2);
+        series_time{i} = local_parse_datetime(string({values.dateTime}'));
+        series_value{i} = str2double({values.value}');
+        keep(i) = true;
+    catch ME
+        warning('MATLAB:read_usgs_file', ...
+            'Failed to process time series %d: %s', i, ME.message);
+    end
 end
 
-times = double(    ...
-     py.mhkit_python_utils.pandas_dataframe.datetime_index_to_ordinal(datapd));
+names = names(keep);
+units = units(keep);
+series_time = series_time(keep);
+series_value = series_value(keep);
+num_series = numel(names);
+
+% Union of all timestamps across parameters, mirroring pandas'
+% combine_first which outer-joins each series on its datetime index
+unique_times = unique(vertcat(series_time{:}));
 
-datast.time = posixtime(datetime(times,                        ...
-                                 'ConvertFrom', 'datenum',     ...
-                                 'TimeZone','UTC'));
+dataset = struct();
+for i = 1:num_series
+    column = nan(numel(unique_times), 1);
+    [is_member, loc] = ismember(series_time{i}, unique_times);
+    column(loc(is_member)) = series_value{i}(is_member);
+
+    dataset.(names(i)) = column;
+    dataset.units.(names(i)) = units(i);
+end
+
+dataset.time = posixtime(unique_times).';
+
+end
+
+function dt = local_parse_datetime(date_strings)
+% Parses USGS dateTime strings (ISO 8601, with or without a UTC offset)
+% into UTC datetimes, mirroring pandas.to_datetime(..., utc=True).
+date_strings = strrep(date_strings, "Z", "+00:00");
+
+has_offset = ~isempty(regexp(date_strings(1), '[+-]\d{2}:\d{2}$', 'once'));
+
+if has_offset
+    dt = datetime(date_strings, ...
+        'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss.SSSX', ...
+        'TimeZone', 'UTC');
+else
+    dt = datetime(date_strings, ...
+        'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss.SSS', ...
+        'TimeZone', 'UTC');
+end
+end
 

mhkit/river/IO/request_usgs_data.m

@@ -37,6 +37,10 @@
 %     write_json : str or None (optional)
 %         Name of json file to write data
 %         to call: request_usgs_data(station,parameter,start_date,end_date,"write_json",write_json)
+%
+%     timeout : double (optional)
+%         Timeout in seconds for each HTTP request. Default = 30
+%         to call: request_usgs_data(station,parameter,start_date,end_date,"timeout",timeout)
 %         
 % Returns
 % -------
@@ -58,10 +62,9 @@
     end_date string
     options.data_type = 'Daily';
     options.write_json string = "";
+    options.timeout (1,1) {mustBeNumeric, mustBePositive} = 30;
 end
 
-MAX_RETRIES = 5;                         % number of query retries if error
-
 % Formulate query
 if strcmp(options.data_type, 'Daily')
     data_url = "https://waterservices.usgs.gov/nwis/dv";
@@ -84,19 +87,8 @@
 disp("Data request URL: " + data_url + api_query)
 
 % Submit query and get data
-for i = 0:MAX_RETRIES
-    try
-        response = webread(data_url + api_query);
-        break;
-    catch ME
-        if i == MAX_RETRIES
-            disp(['MATLAB:request_usgs_data: ', ME.identifier]);
-            rethrow(ME)
-        else
-            pause(1);   % pause(seconds) and retry query
-        end
-    end
-end
+web_options = weboptions('Timeout', options.timeout);
+response = webread_with_retry(data_url + api_query, web_options);
 
 % Parse data label and units
 v = response.value.timeSeries.variable.variableDescription;
@@ -131,3 +123,63 @@
     fid = fopen(options.write_json, 'w');
     fprintf(fid, jsonencode(response));
 end
+
+end
+
+
+function response = webread_with_retry(url, web_options)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Run a USGS web request, retrying transient failures with exponential
+% backoff and jitter
+%
+% Server errors (HTTP 5xx), rate limiting (HTTP 429), request timeouts
+% (HTTP 408), and connection failures are transient, so the request is
+% retried after 1, 2, 4, 8, 16, and 32 seconds, each scaled by a random
+% factor between 0.5 and 1.5 so repeated clients do not retry in
+% lockstep. Other client errors (HTTP 4xx, e.g. an invalid station or
+% parameter) are not transient and are raised immediately.
+%
+% Parameters
+% ------------
+%     url : string
+%         USGS request URL
+%     web_options : weboptions
+%         Options passed to webread
+%
+% Returns
+% ---------
+%     response : struct
+%         Decoded JSON response
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+MAX_RETRIES = 6;
+BASE_DELAY = 1;   % seconds, doubled after each failed attempt
+
+for attempt = 0:MAX_RETRIES
+    try
+        response = webread(url, web_options);
+        return
+    catch ME
+        status = regexp(ME.identifier, 'HTTP(\d+)StatusCodeError', 'tokens', 'once');
+        if isempty(status)
+            is_transient = true;   % timeout or connection failure
+        else
+            status = str2double(status{1});
+            is_transient = status >= 500 || status == 429 || status == 408;
+        end
+
+        if ~is_transient || attempt == MAX_RETRIES
+            disp(['MATLAB:request_usgs_data: ', ME.identifier]);
+            rethrow(ME)
+        end
+
+        delay = BASE_DELAY * 2^attempt * (0.5 + rand());
+        fprintf('USGS request failed (%s), retry %d of %d in %.1f s\n', ...
+            ME.identifier, attempt + 1, MAX_RETRIES, delay);
+        pause(delay);
+    end
+end
+end

mhkit/river/performance/circular.m

@@ -1,29 +1,34 @@
-function [D_E,projected_capture_area]=circular(diameter)
+function [D_E, projected_capture_area] = circular(diameter)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculates the equivalent diameter and projected capture area of a
-%     circular turbine
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Calculates the equivalent diameter and projected capture area of a
+% circular turbine
 %
 % Parameters
 % ------------
-%     diameter : float
-%         Turbine diameter [m]
+% diameter : double [m]
+%   Turbine diameter
 %
 % Returns
 % ---------
-%     D_E : float
-%        Equivalent diameter [m]
-%
-%     projected_capture_area : float
-%         Projected capture area [m^2]
+% D_E : double [m]
+%   Equivalent diameter
+% projected_capture_area : double [m^2]
+%   Projected capture area
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
+arguments (Input)
+    diameter (1,1) {mustBeNumeric}
+end
 
-result=py.mhkit.river.performance.circular(diameter);
+arguments (Output)
+    D_E (1,1) {mustBeNumeric}
+    projected_capture_area (1,1) {mustBeNumeric}
+end
 
-resultc=cell(result);
-D_E=resultc{1};
-projected_capture_area=resultc{2};
+D_E = diameter;
+projected_capture_area = (1/4) * pi * (D_E.^2);
 
+end

mhkit/river/performance/ducted.m

@@ -1,28 +1,34 @@
-function [D_E,projected_capture_area]=ducted(diameter)
+function [D_E, projected_capture_area] = ducted(diameter)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculates the equivalent diameter and projected capture area of a
-%     ducted turbine
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Calculates the equivalent diameter and projected capture area of a
+% ducted turbine
 %
 % Parameters
 % ------------
-%     diameter : float
-%         ducted diameter [m]
+% diameter : double [m]
+%   Duct diameter
 %
 % Returns
 % ---------
-%     D_E : float
-%        Equivalent diameter [m]
+% D_E : double [m]
+%   Equivalent diameter
+% projected_capture_area : double [m^2]
+%   Projected capture area
 %
-%     projected_capture_area : float
-%         Projected capture area [m^2]
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
+arguments (Input)
+    diameter (1,1) {mustBeNumeric}
+end
 
-result=py.mhkit.river.performance.ducted(diameter);
+arguments (Output)
+    D_E (1,1) {mustBeNumeric}
+    projected_capture_area (1,1) {mustBeNumeric}
+end
 
-resultc=cell(result);
-D_E=resultc{1};
-projected_capture_area=resultc{2};
+D_E = diameter;
+projected_capture_area = (1/4) * pi * (D_E.^2);
 
+end