From f6f875310e8a894e5be2909403b66ba105f449e5 Mon Sep 17 00:00:00 2001 From: William Pringle Date: Mon, 16 Mar 2020 22:58:48 -0500 Subject: [PATCH 1/3] allowing for backupdem to replace dem outside of dem bbox If the bbox of geodata is set to be greater than the extents of the main dem, backupdem can replace those values: - geodata now does not set obj.x0y0 to that of the main dem if backupdem present - in edgefx add line to replace nans in main dem with the backupdem (will occur if outside main dem extents) - add line to handle case of fl being a positive scalar to indicate low pass filter. --- @edgefx/edgefx.m | 47 +++++++++++++++++++++++++++++++++++----------- @geodata/geodata.m | 13 ++++++++++--- 2 files changed, 46 insertions(+), 14 deletions(-) diff --git a/@edgefx/edgefx.m b/@edgefx/edgefx.m index 16aabb8f..e9452db5 100644 --- a/@edgefx/edgefx.m +++ b/@edgefx/edgefx.m @@ -399,6 +399,11 @@ % interpolate DEM's bathy linearly onto our edgefunction grid. [xg,yg] = CreateStructGrid(obj); tmpz = feat.Fb(xg,yg); + if ~isempty(feat.Fb2) + tmpz2 = feat.Fb2(xg,yg); + tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz)); + clear tmpz2 + end % Initialise wld obj.wld = NaN([obj.nx,obj.ny]); grav = 9.807; @@ -433,6 +438,11 @@ [xg,yg] = CreateStructGrid(obj); tmpz = feat.Fb(xg,yg); + if ~isempty(feat.Fb2) + tmpz2 = feat.Fb2(xg,yg); + tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz)); + clear tmpz2 + end tmpz(tmpz > 50) = 50; % ensure no larger than 50 m above land % use a harvestine assumption dx = obj.h0*cosd(min(yg(1,:),85)); % for gradient function @@ -465,12 +475,12 @@ else filtit = 0; for lambda = obj.fl' - if all(lambda ~= 0) - % do a bandpass filter - tmpz_ft = filt2(tmpz,dy,lambda,'bp') ; - elseif lambda(2) == 0 + if length(lambda) == 1 || lambda(2) == 0 % do a low pass filter tmpz_ft = filt2(tmpz,dy,lambda(1),'lp') ; + elseif all(lambda ~= 0) + % do a bandpass filter + tmpz_ft = filt2(tmpz,dy,lambda,'bp') ; else % highpass filter not recommended warning(['Highpass filter on bathymetry in slope' ... @@ -636,7 +646,7 @@ tempbb{jj} = feat.boubox; end - [xg,yg]=CreateStructGrid(obj); + [xg,yg] = CreateStructGrid(obj); % STEP 2: For each channel point, set the resolution around a % neighborhood of points as |h|/ch where ch is non-dim # between 0.5-1.5 @@ -670,7 +680,13 @@ end end toc - dp = max(1,-feat.Fb(xg(tidx),yg(tidx))); + tmpz = feat.Fb(xg(tidx),yg(tidx)); + if ~isempty(feat.Fb2) + tmpz2 = feat.Fb2(xg(tidx),yg(tidx)); + tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz)); + clear tmpz2 + end + dp = max(1,-tmpz); obj.chd(tidx) = dp/obj.ch; obj.chd(obj.chd < obj.min_el_ch) = obj.min_el_ch; clearvars Fb pts dp radii tempbb xg yg @@ -804,6 +820,17 @@ function plot(obj,type) hh_m(nearshore & hh_m > obj.max_el_ns) = obj.max_el_ns; end hh_m(hh_m < obj.h0 ) = obj.h0; + + % Compute tmpz if necessary + if length(obj.max_el) > 1 || length(obj.g) > 1 || obj.dt(1) >= 0 + tmpz = feat.Fb(xg,yg); + if ~isempty(feat.Fb2) + tmpz2 = feat.Fb2(xg,yg); + tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz)); + clear tmpz2 + end + end + for param = obj.max_el' if numel(param)==1 && param~=0 mx = obj.max_el(1); @@ -815,8 +842,8 @@ function plot(obj,type) dp1 = param(2); dp2 = param(3); - limidx = (feat.Fb(xg,yg) < dp1 & ... - feat.Fb(xg,yg) > dp2) & (hh_m > mx | isnan(hh_m)); + limidx = (tmpz < dp1 & tmpz > dp2) & ... + (hh_m > mx | isnan(hh_m)); hh_m( limidx ) = mx; end @@ -858,8 +885,7 @@ function plot(obj,type) dp1 = param(2); dp2 = param(3); - limidx = (feat.Fb(xg,yg) < dp1 & ... - feat.Fb(xg,yg) > dp2) ; + limidx = (tmpz < dp1 & tmpz > dp2) ; dmy( limidx ) = lim; end @@ -895,7 +921,6 @@ function plot(obj,type) end % Estimate wavespeed in ocean (second term represents orbital % velocity at 0 degree phase for 1-m amp. wave). - tmpz = feat.Fb(xg,yg); grav = 9.807; % Limit the minimum depth to 1 m (ignore overland) tmpz(tmpz > - 1) = -1; diff --git a/@geodata/geodata.m b/@geodata/geodata.m index 9de8994d..a4ca80be 100644 --- a/@geodata/geodata.m +++ b/@geodata/geodata.m @@ -306,6 +306,8 @@ obj = ParseDEM(obj) ; + + end function obj = ParseShoreline(obj) @@ -556,9 +558,14 @@ clear x y demz else % main interpolant - obj.Fb = griddedInterpolant({x,y},demz,... - 'linear','nearest'); - obj.x0y0 = [x(1),y(1)]; + if obj.BACKUPdemfile~=0 + obj.Fb = griddedInterpolant({x,y},demz,... + 'linear','none'); % no extrapolation (so use Fb2) + else + obj.Fb = griddedInterpolant({x,y},demz,... + 'linear','nearest'); + obj.x0y0 = [x(1),y(1)]; + end % clear data from memory clear x y demz end From b2d23d9db376ecb0150249386ba77099970fd874 Mon Sep 17 00:00:00 2001 From: William Pringle Date: Tue, 17 Mar 2020 09:57:05 -0500 Subject: [PATCH 2/3] updatign msh.interp help for fillinside --- @msh/msh.m | 3 ++- @msh/private/GridData.m | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/@msh/msh.m b/@msh/msh.m index 356eeb1e..1fbb9b71 100644 --- a/@msh/msh.m +++ b/@msh/msh.m @@ -871,7 +871,8 @@ function write(obj,fname,type) % relevant for CA interpolation method). % default value N=1. % - % nan - 'fill' to fill in any NaNs appearing in bathy + % nan - 'fill' to fill in any NaNs everywhere + % -'fillinside' to fill NaNs only in DEM extents % % mindepth - ensure the minimum depth is bounded in the % interpolated region diff --git a/@msh/private/GridData.m b/@msh/private/GridData.m index 0a673f05..21a58b74 100644 --- a/@msh/private/GridData.m +++ b/@msh/private/GridData.m @@ -30,7 +30,8 @@ % relevant for CA interpolation method). % default value N=1. % -% nan (optional) - 'fill' to fill in any NaNs appearing in bathy +% nan (optional) - 'fill' to fill in any NaNs everywhere +% -'fillinside' to fill NaNs only in DEM extents % % mindepth (optional) - ensure the minimum depth is bounded in the % interpolated region From c5b94387c1b68b8ec6b6794543314e867f9aa1fc Mon Sep 17 00:00:00 2001 From: William Pringle Date: Tue, 17 Mar 2020 17:48:54 -0500 Subject: [PATCH 3/3] Update Calc_f13.m adding advection_state --- utilities/Calc_f13.m | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/utilities/Calc_f13.m b/utilities/Calc_f13.m index 5bc01a54..36d658ec 100644 --- a/utilities/Calc_f13.m +++ b/utilities/Calc_f13.m @@ -11,6 +11,7 @@ % 'Mn' ('mannings_n_at_sea_floor') % 'Ss' ('surface_submergence_state') % 'Re' ('initial_river_elevation') +% 'Ad' ('advection_state') % % 3) then either: % 'inpoly' followed by... @@ -47,6 +48,9 @@ elseif strcmpi(attribute,'Re') attrname = 'initial_river_elevation'; default_val = 0; +elseif strcmpi(attribute,'Ad') + attrname = 'advection_state'; + default_val = -999; else error(['Attribute ' attribute ' not currently supported']) end