Adding threshold option to extract_subdomain

README.md

@@ -136,11 +136,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