Shift in Bathymetry when Generating Mesh

[PatMottmac]

Hi All, First time working with this tool and am running into a slight shift in my DEM data during mesh generation. Images below show (1) the raw bathymetry generated by the matlab code and (2) actual DEM data interpolated to the mesh in SMS. Wondering if there's a projection issue in my code that I may be missing? For reference the DEM is in Hor.: NAD83 decimal degrees, Vert.: NAVD88 m.

(1) - Raw bathy output from matlab code. Thalweg input shown as red line:

(2) - DEM Interpolated to mesh during postprocessing. Note actual location of channel:

Matlab code below for reference - Adapted from Example 6

% Example_6_GBAY: Mesh the Galveston bay (GBAY) region in
% high resolution.

clearvars; clc;

addpath('..')
addpath(genpath('../utilities/'))
addpath(genpath('../datasets/'))
addpath(genpath('../m_map/'))

%% STEP 1: set mesh extents and set parameters for mesh.
bbox = [-94.08 -93.65;
    29.5  29.99];
min_el    = 20;         % Minimum mesh resolution in meters.
max_el    = 1000;       % Maximum mesh resolution in meters.
max_el_ns = 100;        % Maximum mesh resolution in nearshore (0.1 deg of shoreline) in meters
min_el_ch = 25;         % Minimum resolution near channel (100 default)
ch        =  0.4;       % How mesh resolution scales with depth
fs        =  4;         % Minimum number of elements to resolve a feature
g         =  0.2;       % Mesh grade in decimal percent - default 0.25.  Higher number equals faster grading from small cells to big
dt        =  2;         %Ensure mesh is stable at this timestep [s]
%Mesh building parametesr
mindepth = 0;           % Minimum depth in mesh
name = 'TexasPB_NEST';  % Mesh name
%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
%coastline = 'us_medium_shoreline_polygon';
coastline = 'TXPB_DEM_0m';
demfile   = 'CUDEM_TXPB.nc';
gdat = geodata('shp',coastline,...
               'dem',demfile,...
               'bbox',bbox,...
               'h0',min_el);
%load ECGC_Thalwegs.mat % Load the Channel thalweg data as pts2 - USER NOTE - can specify your own channels as an XY cell pair
load sabine_thalweg.mat; % Load the Channel thalweg data as pts2 - USER NOTE - can specify your own channels as an XY cell pair
sabine_thalweg{1} = [sabine_thalweg{1} ; [NaN,NaN]];
%% STEP 3: create an edge function class
fh = edgefx('geodata',gdat,...
    'fs',fs,...
    'g',g,...
    'channels',sabine_thalweg,...
    'min_el_ch',min_el_ch,...
    'dt',dt,...
    'max_el',max_el,...
    'ch',ch,...
    'max_el_ns',max_el_ns); 
%% Repeat STEPS 1?3 for a high resolution domain in selected area
min_el = 20; % minimum resolution in meters.
max_el = 40; % maximum resolution in meters.
max_el_ns = 100; % maximum resolution nearshore.
%Bounding box for refined area
bbox = [-94.94 -93.838;
    29.66  29.69];
%Specify coastline and updated dem
coastline = 'us_medium_shoreline_polygon';
demfile   = 'CUDEM_TXPB.nc';
%Create now geographic class
gdat2 = geodata('shp',coastline,'dem',demfile,'h0',min_el);
%Edge function class for refined area
fh2 = edgefx('geodata',gdat2,...
    'fs',fs,...
    'ch',ch,...
    'dt',dt,...
    'channels',sabine_thalweg,...
    'min_el_ch',min_el_ch,...
    'max_el',max_el,...
    'max_el_ns',max_el_ns,...
    'g',g);
%% STEP 4: Pass your edgefx class object along with some meshing options and
% build the mesh...
mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1,'proj','lambert');
% now build the mesh with your options and the edge function.
mshopts = mshopts.build;
%% STEP 5: Plot it and write a triangulation fort.14 compliant file to disk.
m = mshopts.grd;
%m = interp(m,gdat,'mindepth',mindepth,'nan','fill'); % interpolate bathy to the mesh
%m = interp(m,gdat,'mindepth',mindepth,'nan','fillinside',); % interpolate bathy to the mesh
m = interp(m,demfile,'mindepth',mindepth,'nan','fillinside'); % interpolate bathy to the mesh
m = make_bc(m,'auto',gdat); % make the nodestring boundary conditions
% Plotting and writing 
plot(m,'type','bmesh'); % plot bathy on the mesh
plot(m,'type','bdnotri','holdon',1);  % plot the boundary conditions on top of the bathy mesh
write(m,name);

Thanks for the help and for this great tool!

[krober10nd]

Hey @PatMottmac Thanks for the code! Sorry you're experiencing this.

Could you try a different interpolation method (without the fillinside option) and update us what happens?

m = interp(m,demfile,'interp','linear'); 

Either way, you're doing everything right. It looks like a bug, would it be possible to share the DEM (if necessary privately) so I can analyze what's going on?

[PatMottmac]

Thanks @krober10nd ! That seemed to fix it - guess it was just the interpolation. Thanks again for the help

[krober10nd]

Just curious, what version of OM are you using? If you're using version control, you type git status in the folder with it.

I believe someone else encountered this with the CA interp method when interpolating from a DEM with variable resolution horizontal grid spacing (i.e., x grid and y-grid spaces were not equal).

We put in a fix in https://github.com/CHLNDDEV/OceanMesh2D/commit/8311513226f75f6c4e085eeca375324fd1fa897d on Dec 2.

[PatMottmac]

I was using the v4.0 - just downloaded the latest zip file on 3/14/2020, not using version control unfortunately (will switch now -just wanted to get this up and running). My DEM does have variable resolution so maybe that triggered the issue?

[krober10nd]

Ah okay. Yes, it’s likely that we have not fully resolved the non regular DEM grid spacing case.

Would it be possible to have access to the DEM for testing purposes to figure this coding bug out? That would greatly help us figure this out. Thanks

[krober10nd]

Hey Patrick, in #204 there's more information why your DEM didn't work and a potential fix.