Add a tool to remove pits from the mesh bathymetry

[HamishB]

Pits in the bathymetry can make the wetting/drying portion go unstable.
Here's a little code contribution to replace pits with the average bathymetry
of the neighboring nodes, if it may be useful.

% find pits in OceanMesh2D mesh bathymetry
%   Hamish Bowman 10 Sept. 2020

m = msh('fort14_with_pits');
m_orig.b = m.b;

% find node points likely to see some wetting and drying
% may want to get rid of upper bound (here 5m above MSL)
intertidal = find(m.b < 3 & m.b > -5);

sz = size(m.t);
is_pit = zeros(size(m.b));  % init array
mean_neighbors = nan(size(m.b));  % init array

% brute force
% note: this modifies along the way, not atomically at the end
%    so pits may come and go as the neighboring pit nodes get filled
tic
for n = [intertidal']
   idx = find(m.t == n);
   [rows cols] = ind2sub(sz, idx);

   hits = m.t(rows, :);
   neighbors = unique(hits(find(hits ~= n)));

   % ignore spur nodes along the boundary (maybe skip all bdy nodes?)
   %   or run QA code to nuke spur elements first?
   if(length(neighbors) > 2 & m.b(n) > max(m.b(neighbors)))
      disp(['node ' num2str(n) ' (' num2str(m.b(n), '%.3f') ...
         ') is an intertidal pit with ' ...
         num2str(length(neighbors)) ' neighbors who average ' ...
	 num2str(mean(m.b(neighbors)), '%.3f') ' meters depth'])
      is_pit(n) = 1;
      mean_neighbors(n) = mean(m.b(neighbors));

      % replace if pit is more than 10 cm deep (does size of pit matter?)
      if(m.b(n) - max(m.b(neighbors)) > 0.10)
         disp(['  * diff: ' num2str(m.b(n) - max(m.b(neighbors)))])
         m.b(n) = mean(m.b(neighbors));
      end
   end
end
toc

disp([num2str(length(find(is_pit))) ' pits found in mesh bathymetry'])


% or to fill in pits after the scan,  (e.g. if scan is multi-threaded)
%for n = [find(is_pit)']
%   % replace if pit is more than 10 cm deep (does size of pit matter?)
%   if(m.b(n) - mean_neighbors(n) > 0.10)
%      m.b(n) = mean_neighbors(n);
%   end
%end


%%% you may want to iterate through this 2-10 times.

[jasonfleming]

This looks interesting. This seems to be a common issue. Nate Dill has a Perl script for filling depressions also: https://github.com/natedill/ourPerl/blob/master/AdcircUtils/depressionFiller.pl ... I wonder if there are additional codes out there to address this.

[krober10nd]

thanks @HamishB perhaps we should create a folder called "mesh_checks" and put these one off functions there?

[HamishB]

Just to note some tools for filling in pits in raster DEMs, but due to aliasing perhaps it is better to attack the problem after the triangular mesh is created, at least for the case of single node pits.

https://grass.osgeo.org/grass78/manuals/r.fill.dir.html
https://grass.osgeo.org/grass78/manuals/addons/r.hydrodem.html

I see Nate's code goes after peaks as well as depressions. Are small gradient peaks a problem for ADCIRC too?

[jasonfleming]

This is interesting @HamishB ... I'm not an expert on DEMs or meshing, so I don't know the rationale for going after small gradient peaks. I just wanted to chime in with something that would hopefully be helpful :-) (like these additional GRASS links). Seems like depression filling could be its own subsystem in the mesh quality assurance process ...

[krober10nd]

Yea, @jasonfleming we found that depressions sometimes occur when the mesh resolution is at the scale of man-made structures (~1-10 m) like underpasses, or for example a dredged dock that has a thin separation between the two piers. Interpolation can lead to pockets.

Philosophically, I rather edit the mesh than edit the DEM (the source of information).

[krober10nd]

I cleaned this up a little bit and will put it in a PR shortly as a new function of the msh class.

function [m] = fill_bathymetric_holes(m,hole_thres)
% obj = fill_bathymetric_holes(msh_object,hole_thres);
% fill_bathymetric_holes: This function removes holes that 
% appeared perhaps during interpolation of bathymetric data onto the mesh vertices 
% using for example msh.interp (which is a cell-averaging approach). 
%
% A "hole" or "pit" is defined as a node with a neighboring vertex that has a
% bathymetry greater than hole_thres meters and is located in the meshes 
% intertidal zone (+3 to -5 m around datum). This script fills the value 
% at the hole with the average of its neighbors effectively lifting the 
% value up (e.g., Laplacian smoothing). As such, this script can be repeated
% until the desired effect is achieved. 
%
% Input:          m - a mesh object with a bathymetry set defined at its vertices 
%        hole_thres - the differetial in bathymetry between a vertices 
%                     neighobrs that flags the node as a hole 
% Output:         m - a new mesh object with the holes filled in the
%                     bathymetry field of the mesh object.
% 
% Original author: Hamish Bowman Sept. 10, 2020
% Modifications by Keith Roberts Nov. 9, 2020
%
% find node points likely to see some wetting and drying
% may want to get rid of upper bound (here 5m above MSL)
intertidal = find(m.b < 3 & m.b > -5);

sz = size(m.t);
is_pit = zeros(size(m.b));  % init array
mean_neighbors = nan(size(m.b));  % init array

% brute force
% note: this modifies along the way, not atomically at the end
%    so pits may come and go as the neighboring pit nodes get filled
tic
for n = [intertidal']
   idx = find(m.t == n);
   [rows, ~] = ind2sub(sz, idx);

   hits = m.t(rows, :);
   neighbors = unique(hits(find(hits ~= n)));

   % ignore spur nodes along the boundary (maybe skip all bdy nodes?)
   %   or run QA code to nuke spur elements first?
   if(length(neighbors) > 2 & m.b(n) > max(m.b(neighbors)))
      disp(['node ' num2str(n) ' (' num2str(m.b(n), '%.3f') ...
         ') is an intertidal pit with ' ...
         num2str(length(neighbors)) ' neighbors who average ' ...
	 num2str(mean(m.b(neighbors)), '%.3f') ' meters depth'])
      is_pit(n) = 1;
      mean_neighbors(n) = mean(m.b(neighbors));

      % replace if pit is more than hole_thres cm deep (does size of pit matter?)
      if(m.b(n) - max(m.b(neighbors)) > hole_thres)
         disp(['  * diff: ' num2str(m.b(n) - max(m.b(neighbors)))])
         m.b(n) = mean(m.b(neighbors));
      end
   end
end
toc

disp([num2str(length(find(is_pit))) ' pits found in mesh bathymetry'])

[HamishB]

Thanks. Will continue to think about finding & dealing with pits of two or more nodes. A simple way after a tides-only ADCIRC run is to filter the maxele.63 file for nodes where the max elevation is exactly 0.0. (will not catch above MSL depressions, but anecdotally I don't think those are as much trouble as inland ponds with depths below MSL?)

[krober10nd]

Note that using the original DEM resolution by creating a separate geodata instance and using that to do the interpolation of bathymetry seemed to avoid 99% of the holes. But yea, this could still be useful as a checker. Sometimes data has flaws like this and it'd be good to know right away.