Adding threshold option to extract_subdomain (#177)

* Adding depth-based bounds option to extract_subdomain

Author: Keith Roberts

Examples/Example_11_Remeshing_Patches.m

@@ -1,5 +1,5 @@
 % Example_11_Remeshing_Patches.m
-% 
+%
 % This example highlights a workflow on how to re-mesh and insert patches
 % using user defined mesh sizing functions from OceanMesh2D preserving
 % the pathces subdomain boundaries exactly.
@@ -38,7 +38,7 @@
 % Get out the msh class and put on nodestrings
 m = mshopts.grd;
 % Interpolate topobathy data onto the vertices of the mesh.
-m = interp(m,'SRTM15+V2.1.nc'); 
+m = interp(m,'SRTM15+V2.1.nc');
 %% Extract the region which you want to remesh
 % For the purpose of this example, we have drawn a random polygon on top
 % of the mesh for which we would like to remesh.
@@ -50,27 +50,27 @@
   173.5714  -43.6053
   173.0915  -44.1310];
 % This extracts the parent mesh but with the polygon removed.
-subdomain = extract_subdomain(m,hole); 
+subdomain = extract_subdomain(m,hole);
 % Visualize the subdomain
 plot(subdomain);
 %% Get the boundary of this hole in the mesh
-% Follow the instructions in the title of the plot and click on one of the 
+% Follow the instructions in the title of the plot and click on one of the
 % points on the hole's boundary and type the number into the screen.
 poly = get_poly(subdomain);
-%% Re-meshing 
-% Using the polygon we just extracted (poly) and the original mesh size 
-% function (fh), remesh this hole with Delaunay refinement mesh2d. 
+%% Re-meshing
+% Using the polygon we just extracted (poly) and the original mesh size
+% function (fh), remesh this hole with Delaunay refinement mesh2d.
 % Recall poly is a cell-array so pass it the hole you want to mesh!
-subdomain_new = mesh2dgen(poly{1},fh); 
+subdomain_new = mesh2dgen(poly{1},fh);
 %% Merge the patch in
 % Remove the subdomain from the parent mesh
-m_w_hole = extract_subdomain(m,hole,1);
+m_w_hole = extract_subdomain(m,hole,"keep_inverse",1);
 % Since the boundaries match identically, this is a trivial merge.
 m_new = plus(subdomain_new, m_w_hole, 'match');
-% Use a trivial nearest neighbor interpolation to put the bathy back on 
+% Use a trivial nearest neighbor interpolation to put the bathy back on
 % the new mesh
 ind = nearest_neighbor_map(m, m_new);
-m_new.b = m.b(ind); 
+m_new.b = m.b(ind);
 % Plot the bathy on the mesh with hole polygon in red
-plot(m_new,'type','bmesh'); hold on; 
+plot(m_new,'type','bmesh'); hold on;
 plot(poly{1}(1:end-1,1),poly{1}(1:end-1,2),'r-','linewi',3);

README.md

@@ -138,11 +138,12 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 ### Unreleased
 
 ## Added
-- `mesh2d` interface improvements to filter small polygons. 
+- `mesh2d` interface improvements to filter small polygons.
 - Support for creation of `fort.20` files for forcing rivers by @Jiangchao3
 
 ## Changed
 - `msh.plot()` overhaul. All options specified via kwarg.
+- `utilities/extract_subdomain` now is called with kwargs.
 
 ## Fixed
 - Boundary labeling fix

utilities/extract_subdomain.m

@@ -1,45 +1,63 @@
-function [obj,ind] = extract_subdomain(obj,bou,keep_inverse,centroid,nscreen)
-% [obj,ind] = extract_subdomain(obj,bou,keep_inverse,centroid,nscreen)
-% 
+function [obj,ind] = extract_subdomain(obj,bou,varargin)
+% [obj,ind] = extract_subdomain(obj,bou)
+%
 % Inputs:
 % bou: a bbox, i.e.: [lon min, lon_max;
-%                     lat_min, lat_max], or    
+%                     lat_min, lat_max], or
 %      a NaN-delimited polygon of the domain to extract
 % keep_inverse: = 0 [default] to get the sub-domain inside the bou polygon
 %               = 1 to get the sub-domain outside the bou polygon
 % centroid: = 0 [default] inpolygon test is based on whether all vertices
 %              of the element are inside (outside) the bou polygon
 %           = 1 inpolygon test is based on whether the element centroid
 %              is inside (outside) the bou polygon
+% min_depth = topobathymetry value to keep elements inside bou so that they have a
+%             sufficiently deep value at their centroid
+% max_depth = topobathymetry value to keep elements inside bou so that they have a
+%             sufficiently shallow value at their centroid.
 % nscreen: = 1 [default] display the notice to screen
 %          = 0 do not display the notice to screen
-% 
+%
 % Outputs:
 % obj: the subset mesh obj (only p and t, properties untouched)
 % ind: an array of indices that can be used to map the mesh properties to
-% the output mesh subset with a subsequent call to "map_mesh_properties". 
+% the output mesh subset with a subsequent call to "map_mesh_properties".
 %
-p = obj.p; t = obj.t; 
-if nargin == 1 || (nargin == 3 && isempty(bou))
+keep_inverse = 0 ;
+centroid = 0 ;
+min_depth = -99999;
+max_depth = +99999;
+nscreen = 1;
+% Otherwise, name value pairs specified.
+% Parse other varargin
+for kk = 1:2:length(varargin)
+    if strcmp(varargin{kk},'keep_inverse')
+        keep_inverse = varargin{kk+1};
+    elseif strcmp(varargin{kk},'centroid')
+        centroid = varargin{kk+1};
+    elseif strcmp(varargin{kk},'nscreen')
+        nscreen = varargin{kk+1};
+    elseif strcmp(varargin{kk},'min_depth')
+        min_depth = varargin{kk+1};
+    elseif strcmp(varargin{kk},'max_depth')
+        max_depth = varargin{kk+1};
+
+    end
+end
+
+p = obj.p; t = obj.t; b = obj.b;
+if nargin == 1 || (nargin == 2 && isempty(bou))
     plot(p(:,1),p(:,2),'k.');
     h = impoly;
     bou  = h.getPosition;
 end
-if nargin < 3 || isempty(keep_inverse)
-    keep_inverse = 0;
-end
-if nargin < 4 || isempty(centroid)
-    centroid = 0;
-end
-if nargin < 5 || isempty(nscreen)
-    nscreen = 1;
-end
+
 % converting bbox to polygon
 if size(bou,1) == 2
      bou = [bou(1,1) bou(2,1);
-            bou(1,1) bou(2,2); 
+            bou(1,1) bou(2,2);
             bou(1,2) bou(2,2);
-            bou(1,2) bou(2,1); 
+            bou(1,2) bou(2,1);
             bou(1,1) bou(2,1)];
 end
 nans = false;
@@ -56,18 +74,23 @@
         in = inpoly(bxyc,bou,edges);
     end
 else
-    bxy1 = p(t(:,1),:); bxy2 = p(t(:,2),:); bxy3 = p(t(:,3),:); 
+    bxy1 = p(t(:,1),:); bxy2 = p(t(:,2),:); bxy3 = p(t(:,3),:);
     if ~nans
-        in1 = inpoly(bxy1,bou); 
-        in2 = inpoly(bxy2,bou); 
+        in1 = inpoly(bxy1,bou);
+        in2 = inpoly(bxy2,bou);
         in3 = inpoly(bxy3,bou);
     else
-        in1 = inpoly(bxy1,bou,edges); 
-        in2 = inpoly(bxy2,bou,edges); 
+        in1 = inpoly(bxy1,bou,edges);
+        in2 = inpoly(bxy2,bou,edges);
         in3 = inpoly(bxy3,bou,edges);
     end
     in = in1 & in2 & in3;
 end
+if min_depth ~= -99999 | max_depth ~= +99999
+     bem = max(b(t),[],2);   % only trim when all vertices
+     selected = bem > min_depth & bem < max_depth;
+     in = logical(in .* selected);
+end
 if keep_inverse == 0
     t(~in,:) = [];
 else
@@ -76,7 +99,7 @@
 % Remove uncessary vertices and reorder triangulation
 [p1,t,ind] = fixmesh(p,t);
 % Put back into the msh obj
-obj.p = p1; obj.t = t;  
+obj.p = p1; obj.t = t;
 %
 if nscreen
     % Displaying notice for mapping mesh properties