Expanding Manning's landcover option (CCAP and NLCD)

@msh/msh.m

@@ -887,13 +887,13 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
             %                         default value N=1.
             %
             %        nan - 'fill' to fill in any NaNs everywhere
-            %            -'fillinside' to fill NaNs only in DEM extents
+            %            - 'fillinside' to fill NaNs only in DEM extents
             %
             %   mindepth - ensure the minimum depth is bounded in the
-            %                         interpolated region
+            %              interpolated region
             %
             %   maxdepth - ensure the maximum depth is bounded in the
-            %                         interpolated region
+            %              interpolated region
             %
             %   ignoreOL - = 0 [default]: interpolate data without masking overland
             %              = 1          : Mask overland data which may help for seabed-only interpolation
@@ -903,7 +903,9 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
             %                             (i.e., keeps the sign of the topography the same as the dem)
             %              = 1 [default]: inverts the values of the dem
             %                             (underwater is positive depth)
-
+            %   lut      - A look up table (lut). See nlcd and ccap in
+            %              datasets/ for examples
+
             % if give cell of geodata or dems then interpolate all
             if iscell(geodata) || isstring(geodata)
                 for i = 1:length(geodata)

@msh/private/GridData.m

@@ -40,6 +40,9 @@
 %                         interpolated region 
 %
 %   ignoreOL (optional) - NaN overland data for more accurate seabed interpolation
+%   
+%   lut (optional)      - A look up table (lut). See nlcd and ccap in
+%                         datasets/ for examples
 %
 % slope_calc (optional) - 'rms' [default]: Compute cell-averaged slope based on root-mean-square
 %                       - 'abs'          : Compute cell-averaged slope based on mean of absolute value
@@ -62,7 +65,7 @@
 %   regions
 %
 %   Edits by Keith Roberts, 2019-4-4 to interpolate the floodplain
-
+%
 %% Test optional arguments
 % default
 nn = length(obj.p);
@@ -76,9 +79,10 @@
 maxdepth = +inf ; 
 rms_slope_calc = true;
 slope_calc = 'rms';
+lut = [] ; 
 if ~isempty(varargin)
     varargin=varargin{1} ; 
-    names = {'K','type','interp','nan','N','mindepth','maxdepth','ignoreOL','slope_calc'};
+    names = {'K','type','interp','nan','N','mindepth','maxdepth','ignoreOL','slope_calc','lut'};
     for ii = 1:length(names)
         ind = find(~cellfun(@isempty,strfind(varargin(1:2:end),names{ii})));
         if ~isempty(ind)
@@ -111,6 +115,8 @@
                 else
                    error(['Invalid entry, ' temp ', for slope_calc'])
                 end
+            elseif ii == 10
+                lut = varargin{ind*2};
             end
         end    
     end
@@ -150,6 +156,10 @@
             'set. Change nan to "fillinside" to avoid this.'])
 end
 
+if ~isempty(lut)
+    disp('Using look up table...'); 
+end
+
 %% Let's read the LON LAT of DEM if not already geodata
 flipUD = 0;
 if ~isa(geodata,'geodata')
@@ -302,6 +312,12 @@
 % bound all depths above maxdepth 
 DEM_Z(DEM_Z > maxdepth) = maxdepth ; 
 
+% if look up table is passed
+if ~isempty(lut)
+    DEM_Z(isnan(DEM_Z) | DEM_Z == 0)=length(lut); % <--default value to NaN
+    % Convert to Mannings
+    DEM_Z = lut(abs(round(DEM_Z)));
+end
 %% Make the new bx, by, b and calculate gradients if necessary
 if strcmp(type,'slope') || strcmp(type,'all')
     bx = NaN(length(K),1); 

Examples/Example_6b_GBAY_w_floodplain.m

@@ -50,6 +50,7 @@
 % 10-m contour extracted from the Coastal Relief Model.
 coastline = 'us_coastalreliefmodel_10mLMSL';
 demfile   = 'galveston_13_mhw_2007.nc';
+landuse   = 'galveston_2016_CCAP.nc'; 
 
 gdat = geodata('shp',coastline,...
                'dem',demfile,...
@@ -69,7 +70,12 @@
 % plot resolution on the mesh
 plot(m,'type','resomesh','colormap',[10 0 1e3])
 
-%interpolate bathy using special constraining technique for overland and
-%underwater
+% interpolate bathy using special constraining technique
+% for overland and underwater
 m = interpFP(m,gdat,muw,gdatuw);
 plot(m,'type','bmesh') % plot the bathy on the mesh
+
+% computing mannings based on CCAP landcover data using 
+% cell-averaged interpolation with stencil 4*grid_size for stability
+m = Calc_Mannings_Landcover(m,landuse,'ccap','N',4);
+plot(m,'type','mann') % plot the mannings on the mesh

README.md

@@ -154,6 +154,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 - Correctly deleting weirs from boundary object through `make_bc` delete method. See https://github.com/CHLNDDEV/OceanMesh2D/pull/205
 - Array format fix for reading in ibtype and nvell from fort.14 file and when executing carry_over_weirs. See https://github.com/CHLNDDEV/OceanMesh2D/pull/206
 - Fix for irregular grid spacings in DEMs. See https://github.com/CHLNDDEV/OceanMesh2D/pull/204
+- tidal constituents for `Make_f15` can now contain "major8" in addition to other constituents in the string/cell array https://github.com/CHLNDDEV/OceanMesh2D/pull/221
+## Changed
+- Renamed `Calc_NLCD_Mannings` to `Calc_Mannings_Landcover` and making option for 'ccap' landcover type in addition to 'nlcd' (default) and added the ability to using user specified inteprolation (e.g., nearest, linear, cell-averaging, etc.) of the landcover data to the mesh vertices. https://github.com/CHLNDDEV/OceanMesh2D/pull/221
 
 ### [4.0.0] - 2021-03-14
 ## Added

utilities/Calc_Mannings_Landcover.m

@@ -0,0 +1,93 @@
+function obj = Calc_Mannings_Landcover(obj,data,type,varargin)
+% obj = Calc_Mannings_Landcover(obj,data,type,varargin)
+% Input a msh class object and interpolates a land-cover database file
+% (netcdf format) onto the msh while converting to mannings values through
+% a look-up table
+% 
+% type:
+%  - 'nlcd': 
+%    NLCD 2006 Land Cover (2011 Edition) (1.1Gb) from
+%    https://www.mrlc.gov/nlcd06_data.php
+%  - 'ccap':
+%    CCAP NOAA land cover:
+%    https://coast.noaa.gov/data/digitalcoast/pdf/ccap-class-scheme-regional.pdf
+%
+% varargin: Accepts the same options as for msh.interp to control how 
+%           data is interpolated; see 'help msh.interp'
+% 
+%  Author:            WP July 19, 2018
+%  Update for CCAP:   WP May 5, 2021
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if nargin < 3 || isempty(type)
+   type = 'nlcd'; 
+end
+
+if strcmp(type,'nlcd')
+    disp('Info: Using NLCD table')
+    % NLCD table
+    load nlcd
+    nlcd_mannings(end+1)=nlcd_mannings(11); % <--make NaN = default value
+    varargin{end+1}='lut';
+    varargin{end+1}=nlcd_mannings;
+elseif strcmp(type,'ccap')
+    disp('Info: Using CCAP table')
+    % CCAP table
+    load ccap 
+    ccap_mannings(end+1)=ccap_mannings(21); % <-- make NaN = default value
+    varargin{end+1}='lut';
+    varargin{end+1}=ccap_mannings;
+else
+    error('Land-cover database not supported')
+end
+% The mannings name and default value
+attrname = 'mannings_n_at_sea_floor';
+default_val = varargin{end}(end); %end of the lut is default
+dmy = msh();  dmy.p = obj.p; dmy.t = obj.t; 
+% Convert to Mannings and interpolate how the user wants
+dmy = interp(dmy,data,varargin{:});
+Man = dmy.b';
+Man( isnan(Man) | Man==0 ) = default_val; %points outside of lut can still be NaN
+%% Make into f13 struct
+if isempty(obj.f13)
+    % Add add mannings as first entry in f13 struct
+    obj.f13.AGRID = obj.title;
+    obj.f13.NumOfNodes = length(obj.p);
+    obj.f13.nAttr = 1;
+    NA = 1;
+else
+    broken = 0;
+    for NA = 1:obj.f13.nAttr
+        if strcmp(attrname,obj.f13.defval.Atr(NA).AttrName)
+            broken = 1;
+            % overwrite existing tau0
+            break
+        end
+    end
+    if ~broken
+        % add mannings to list
+        obj.f13.nAttr = obj.f13.nAttr + 1;
+        NA = obj.f13.nAttr;
+    end
+end
+
+% Default Values
+obj.f13.defval.Atr(NA).AttrName = attrname;
+% We can just put in the options here
+obj.f13.defval.Atr(NA).Unit = 'unitless';
+valpernode = 1;
+obj.f13.defval.Atr(NA).ValuesPerNode = valpernode ;
+obj.f13.defval.Atr(NA).Val = default_val ;
+
+% User Values
+obj.f13.userval.Atr(NA).AttrName = attrname ;
+numnodes = length(find(Man ~= default_val));
+obj.f13.userval.Atr(NA).usernumnodes = numnodes ;
+% Print out list of nodes for each
+K = find(Man ~= default_val);
+obj.f13.userval.Atr(NA).Val = [K ; Man(K)];
+%EOF
+end
+
+
+

utilities/tide_fac.m

@@ -21,8 +21,9 @@
 lat = mean(obj.p(:,2));
 
 %% Allow for an entry just to be major8 to get the major eight constituents
-if strcmp(incnstit{1},'major8')
-   incnstit = {'M2','S2','N2','K2','K1','O1','Q1','P1'};
+jj = strcmp(incnstit,'major8');
+if sum(jj) > 0
+   incnstit = [{'M2','S2','N2','K2','K1','O1','Q1','P1'} incnstit(~jj)];
    disp('Using the major eight harmonic tidal constituents')
    disp(incnstit)
 end
@@ -233,4 +234,4 @@
 ader=[sc;hc;pc;npc;ppc]*dargs./360.0;
 dtau=1.0+ader(2,:)-ader(1,:);
 ader=[dtau;ader];
-end
+end