Closed
Description
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!
Metadata
Metadata
Assignees
Labels
No labels