Skip to content

Commit

Permalink
Merge pull request #61 from CHLNDDEV/dev-testing
Browse files Browse the repository at this point in the history 
Allows for a back-up DEM to replace topo-bathymetric data outside of the main DEM's extent (in the case that the main DEM extent's are not sufficient to cover the entire domain).
  • Loading branch information
keith roberts authored Mar 18, 2020
2 parents 79727f7 + c5b9438 commit b962fc5
Show file tree
Hide file tree
Showing 5 changed files with 54 additions and 16 deletions.
47 changes: 36 additions & 11 deletions @edgefx/edgefx.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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' ...
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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);
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down
13 changes: 10 additions & 3 deletions @geodata/geodata.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -306,6 +306,8 @@

obj = ParseDEM(obj) ;



end

function obj = ParseShoreline(obj)
Expand Down Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion @msh/msh.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion @msh/private/GridData.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions utilities/Calc_f13.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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...
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit b962fc5

Please sign in to comment.