Add additional dolfyn functionality

mhkit/dolfyn/adp/fillgaps_depth.m

@@ -0,0 +1,68 @@
+function out = fillgaps_depth(var, options)
+% Fill gaps (NaN values) in var along depth profile using the specified
+% method in options.method
+% 
+% Parameters
+% ----------
+% var : struct (from a fieldname in dataset in create_dataset function)
+%   The variable to clean
+% method : char
+%   Interpolation method to use. Default is 'spline'
+%   see fillmissing method for other options 
+%   https://www.mathworks.com/help/matlab/ref/fillmissing.html#bva1z1c-method
+% options.maxgap : double 
+%   Maximum gap of missing data to interpolate across. Default is None
+%
+% Returns
+% -------
+% out : struct (from a fieldname in dataset in create_dataset function)
+%  The input struct 'var' with gaps in 'var' interpolated across depth
+%
+% See Also
+% --------
+% create_dataset.m function
+
+% set default value for options.method to 'spline'
+    arguments
+        var;
+        options.method char = 'spline';
+        options.maxgap double = NaN;
+    end
+    % Make sure that var is a structure
+    if ~isstruct(var)
+        ME = MException('MATLAB:fillgaps_depth','var must be a structure');
+        throw(ME);
+    end
+
+    % Make sure that var contains the dolyn fields
+    if ~isfield(var,'coords') || ~isfield(var,'data') 
+        ME = MException('MATLAB:fillgaps_depth',['The provided data ' ...
+            'structure does not appear to have been created by dolfyn']);
+        throw(ME);
+    end
+
+    range_dim = 0;
+    % find first range dim in var.dims
+    for k=1:numel(var.dims) 
+       if contains(var.dims{k}, 'range')
+           % as soon as range_dim is set once, break out of loop
+           range_dim = k;
+           break
+       end
+    end
+
+    if range_dim == 0
+        ME = MException('MATLAB:fillgaps_depth',['No range dimension found']);
+        throw(ME);
+    end
+
+    % now interpolate the nans away, note fillmissing will look for NaNs in
+    % var.data as the missing entries
+    out = var;
+    if isnan(options.maxgap) %if maxgap is Nan, no maxgap was set:
+        out.data = fillmissing(var.data,options.method,range_dim);
+    end
+    if ~isnan(options.maxgap) % if maxgap is not NaN, a maxgap is set:
+        out.data = fillmissing(var.data,options.method,range_dim,MaxGap=options.maxgap);
+    end
+    end

mhkit/dolfyn/adp/fillgaps_time.m

@@ -0,0 +1,68 @@
+function out = fillgaps_time(var, options)
+% Fill gaps (NaN values) in var across time using the specified method in
+% options.method
+% 
+% Parameters
+% ----------
+% var : struct (from a fieldname in dataset in create_dataset function)
+%   The variable to clean
+% options.method : char
+%   Interpolation method to use. Default is 'spline'
+%   see fillmissing method for other options 
+%   https://www.mathworks.com/help/matlab/ref/fillmissing.html#bva1z1c-method
+% options.maxgap : double 
+%   Maximum gap of missing data to interpolate across. Default is None
+%
+% Returns
+% -------
+% out : struct (from a fieldname in dataset in create_dataset function)
+%  The input struct 'var' with gaps in 'var' interpolated across time
+%
+% See Also
+% --------
+% create_dataset.m function
+
+% set default value for options.method to 'spline'
+    arguments
+        var;
+        options.method char = 'spline';
+        options.maxgap double = NaN;
+    end
+    % Make sure that var is a structure
+    if ~isstruct(var)
+        ME = MException('MATLAB:fillgaps_time','var must be a structure');
+        throw(ME);
+    end
+
+    % Make sure that var contains the dolyn fields
+    if ~isfield(var,'coords') || ~isfield(var,'data') 
+        ME = MException('MATLAB:fillgaps_time',['The provided data ' ...
+            'structure does not appear to have been created by dolfyn']);
+        throw(ME);
+    end
+
+    time_dim = 0;
+    % find first time dim in var.dims
+    for k=1:numel(var.dims) 
+       if contains(var.dims{k}, 'time')
+           % as soon as time_dim is set once, break out of loop
+           time_dim = k;
+           break
+       end
+    end
+
+    % if time_dim is still 0, then no time_dim was found in var.dims
+    if time_dim == 0
+        ME = MException('MATLAB:fillgaps_time',['No time dimension found']);
+        throw(ME);
+    end
+
+    % now interpolate the nans away
+    out = var;
+    if isnan(options.maxgap) %if maxgap is Nan, no maxgap was set:
+        out.data = fillmissing(var.data,options.method,time_dim);
+    end
+    if ~isnan(options.maxgap) % if maxgap is not NaN, a maxgap is set:
+        out.data = fillmissing(var.data,options.method,time_dim,'MaxGap', options.maxgap);
+    end
+end

mhkit/dolfyn/adp/find_surface.m

@@ -0,0 +1,83 @@
+function out = find_surface(ds, options)
+%TODO: update comments here since they are just the python comments
+%unchanged
+% Find the surface (water level or seafloor) from amplitude data and
+% adds the variable "depth" to the input Dataset.
+%
+% Parameters
+% ----------
+% vds : struct (from a fieldname in dataset in create_dataset function)
+%   The full adcp dataset
+% options.thresh : int
+%   Specifies the threshold used in detecting the surface. Default = 10
+%     (The amount that amplitude must increase by near the surface for it to
+%     be considered a surface hit)
+% options.nfilt : int
+%   Specifies the width of the median filter applied, must be odd.
+%     Default is None
+%
+% Returns
+% -------
+% None, operates "in place"
+
+    % set default arguments
+    arguments
+        ds;
+        options.thresh {mustBeInteger} = 10;
+        options.nfilt {mustBeInteger} = 0;
+        %TODO check nfilt is odd
+    end
+
+    %TODO/NOTES: So far, this is the best attempt at translating this
+    %function from python into matlab...Since python arrays are 0 indexed
+    %and matlab arrays are 1 indexed some of the indexing in the
+    %translation gets a bit messy. A lot of those indexes need to be double
+    %checked, and there are likely a lot of off-by-one type errors.
+    %Further, this has en error in it at "d1 = ds.coords.range(inds)" that
+    %makes me think some of the indecies need to change even more than 1?
+    %In matlab, the data has dimension Lx1x28xM and I think the second
+    %dimension, 1, gets squeezed out of the python data. In other words,
+    %the line below should read "[~, inds] = max(ds.amp.data, [], 3,
+    %"linear")" to skip over the 1 index. (Not sure on this, but just my
+    %initial thoughts on the first error in this function)
+
+    % This finds the indices of the maximum of the echo profile:
+    [~, inds] = max(ds.amp.data, [], 2, "linear")
+    % This finds the first point that increases (away from the profiler) in
+    % the echo profile
+    % almost the same as python call, just use axis 2 since matlab 1
+    % indexed.
+    edf = diff(cast(ds.amp.data,"int16"), 2);
+    endint = size(ds.vel.data,3);
+    % use matlab : operator to get close to np.arange
+    inds2 = max(uint8(edf<0) .* reshape(uint8(1:endint), 1, [], 1), [], 2);
+    % Calculate the depth of these quantities
+    d1 = ds.coords.range(inds);
+    d2 = ds.coords.range(inds2);
+    % Combine them:
+    D = vertcat(d1, d2);
+    % Take the median value as the estimate of the surface:
+    d = median(D, 1);
+
+    % Throw out values that do not increase near the surface by *thresh*
+    for ip = 1:ds.vel.dims(2)
+        itmp = min(inds(:, ip));
+        if all(edf(itmp:end, :, ip) < thresh)
+            d(ip) = NaN;
+        end
+    end
+
+    if exists('nfilt', 'var') && nfilt ~= 0
+        dfilt = medfilt2(d, [1 nfilt]);
+        dfilt(conv2(isnan(d), ones(1, nfilt) / nfilt, 'same') > 4) = NaN;
+        dfilt(dfilt == 0) = NaN;
+        d = dfilt;
+    end
+
+    ds.depth.data = d;
+    ds.depth.dims = {'time'};
+    ds.depth.units='m';
+
+    out = ds;
+end
+

mhkit/dolfyn/adp/nan_beyond_surface.m

@@ -6,8 +6,8 @@
 % ----------
 % ds : Dataset
 %   The adcp dataset to clean
-% val : nan or numeric
-%   Specifies the value to set the bad values to (default np.nan).
+% val : NaN or numeric
+%   Specifies the value to set the bad values to (default NaN).
 % 
 % Returns
 % -------
@@ -22,7 +22,7 @@
 
     arguments
         ds         
-        options.val = nan;
+        options.val = NaN;
     end
 
     fn = fieldnames(ds);
@@ -56,7 +56,7 @@
         range_limit = ((ds.depth.data-ds.attrs.h_deploy).*cos(beam_angle) ...
                         - ds.attrs.cell_size) + ds.attrs.h_deploy;
     else
-        range_limit = ds.depth.data .* np.cos(beam_angle) - ...
+        range_limit = ds.depth.data .* cos(beam_angle) - ...
             ds.attrs.cell_size;
     end
 

mhkit/dolfyn/adp/val_exceeds_thresh.m

@@ -0,0 +1,48 @@
+function out = val_exceeds_thresh(var, options)
+% Find values of a variable that exceed a threshold value,
+% and assign `options.val` to the velocity data where the threshold is
+% exceeded.
+%
+% Parameters
+% ----------
+% var : struct (from a fieldname in dataset in create_dataset function)
+%   The variable to clean
+% options.thresh : double
+%   The maximum value of data to screen. Default = 5
+% options.val : NaN or double
+%   Specifies the value to set the bad values to. Default is `NaN`
+%
+% Returns
+% -------
+% out : struct (from a fieldname in dataset in create_dataset function)
+%   The input struct `var` with values beyond `options.thresh` set to
+%   `options.val`
+
+    % set default arguments
+    arguments
+        var;
+        options.thresh double = 5.0;
+        options.val double = NaN;
+    end
+
+    % Make sure that var is a structure
+    if ~isstruct(var)
+        ME = MException('MATLAB:val_exceeds_thresh','var must be a structure');
+        throw(ME);
+    end
+
+    % Make sure that var contains the dolyn fields
+    if ~isfield(var,'data')
+        ME = MException('MATLAB:val_exceeds_thresh',['The provided data ' ...
+            'structure does not appear to have been created by dolfyn']);
+        throw(ME);
+    end
+
+    % copy var.data for output
+    out = var;
+    % make logical matrix for if elmnt is beyond options.thresh
+    flags = abs(out.data) > options.thresh;
+    % replace out.data from flags with options.val
+    out.data(flags) = options.val;
+end
+

mhkit/dolfyn/io/read_netcdf.m

@@ -109,8 +109,8 @@
                 else
                     ds.(name).dims{kk} = dimensions(kk).Name;
                 end
-                ds.(name).coords.(dimensions(kk).Name) = ...
-                    ds.coords.(dimensions(kk).Name);
+                ds.(name).coords.(ds.(name).dims{kk}) = ...
+                    ds.coords.(ds.(name).dims{kk});
             end 
             if ~isempty(attrs)
                 for ii = 1:numel(attrs)