incremental Mannings n creation with cell-averaging

[krober10nd]

Describe the bug
When creating the nodal attribute Mannings n landcover parameter, we often have to partition the domain into chunks/slices to fit into available RAM. When this happens, sometimes N>1 cell-averaging stencil overlaps partitions and values can become duplicated leading to errors in the fort.13 file creation for a model such as ADCIRC.

Onur could you please paste here the latest mannings n function we were working on so no one forgets for the next time?

[Onur-Kurum]

Here you go Keith. Please compare with yours. I think this one is not the final working version.

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));
% Print out list of nodes for each
K = find(Man ~= default_val);
if isfield(obj.f13.userval.Atr(NA),'Val') &&  ~isempty(obj.f13.userval.Atr(NA).usernumnodes)
    obj.f13.userval.Atr(NA).usernumnodes = obj.f13.userval.Atr(NA).usernumnodes  + numnodes ;
    obj.f13.userval.Atr(NA).Val = [obj.f13.userval.Atr(NA).Val'; K' , Man(K)']';
else
    obj.f13.userval.Atr(NA).usernumnodes = numnodes ;
    obj.f13.userval.Atr(NA).Val = [K ; Man(K)];
end
%EOF
end

[krober10nd]

Thanks! we'd also like to support custom default values for this calculation.

[krober10nd]

FYI if you have overlapping tiles, the above approach cannot be used.