Merge pull request #204 from CHLNDDEV/interp_irregular_grid_fix

Fix for irregular grid in FindLinearIdx

README.md

@@ -146,6 +146,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 ## Fixed
 - Correctly deleting weirs from boundary object through `make_bc` delete method. See https://github.com/CHLNDDEV/OceanMesh2D/pull/205
 - Array format fix for reading in ibtype and nvell from fort.14 file and when executing carry_over_weirs. See https://github.com/CHLNDDEV/OceanMesh2D/pull/206
+- Fix for irregular grid spacings in DEMs. See https://github.com/CHLNDDEV/OceanMesh2D/pull/204
 
 ### [4.0.0] - 2021-03-14
 ## Added

Tests/RunTests.m

@@ -15,6 +15,8 @@
    TestSanity 
 
    TestEleSizes
+
+   TestInterp
 
    if exist('SRTM15+V2.1.nc')
 

Tests/TestInterp.m

@@ -0,0 +1,71 @@
+%% Test for the msh.interp() method using topo-bathy
+%%  DEMs with following grid spacings: 
+% i)   constant dx = constant dy 
+% ii)  constant dx ~= constant dy 
+% iii) varying dx ~= varying dy
+cd ..
+
+addpath(genpath('utilities/'))
+addpath(genpath('datasets/'))
+addpath(genpath('m_map/'))
+
+% add the geospatial data filenames
+coastline = 'GSHHS_f_L1';
+DEMS = ["CUDEM_equal.nc" "CUDEM_unequal.nc" "CUDEM_varying.nc"];
+
+% meshing parameters
+min_el    = 40;          % minimum resolution in meters.
+max_el    = 500;         % maximum resolution in meters. 
+grade     = 0.25;        % mesh grade in decimal percent. 
+dis       = grade; 	 % increasing resolution with distance from shoreline.
+
+volumes = [];
+% loop over the different DEMs
+for dem = DEMS
+   dem_name = dem{1}; 
+   %x = ncread(dem_name,'lon');
+   %y = ncread(dem_name,'lat');
+   %min(diff(x))
+   %max(diff(x))
+   %min(diff(y))
+   %max(diff(y))
+   %continue  
+
+   gdat = geodata('shp',coastline,'dem',dem_name,'h0',min_el);
+
+   fh = edgefx('geodata',gdat,...
+           'dis',dis,'max_el',max_el,'g',grade);
+
+   mshopts = meshgen('ef',fh,'bou',gdat,...
+                  'plot_on',0,'proj','utm');
+   mshopts = mshopts.build; 
+
+   m = mshopts.grd; 
+
+   m = interp(m,dem_name); 
+
+   %plot(m,'type','bmesh')
+   %print([dem_name(1:end-3) '_bmesh.png'],'-r300','-dpng')
+
+   % compute the volume of the mesh
+   x = m.p(:,1); y = m.p(:,2);
+   areas = polyarea(x(m.t)',y(m.t)');
+   [~,bc] = baryc(m);
+   volumes = [volumes areas*bc];
+end
+
+% Check the different mesh volumes
+disp('Mesh volumes:')
+disp(volumes)
+
+volume_diff = (max(volumes)-min(volumes))/mean(volumes);
+disp('Maximum relative volume difference:')
+disp(volume_diff)
+
+if volume_diff > 0.001
+    error('Mesh volumes are too disparate with different dems')
+    disp('Not Passed: Interp');
+else
+    disp('Mesh volumes are within 0.1% of each other')
+    disp('Passed: Interp');
+end

utilities/FindLinearIdx.m

@@ -4,32 +4,55 @@
 % lon and lat must be matrices created by ndgrid. 
 % kjr, 20171210 chl, und.
 % wjp, 20201120 making sure dx and dy can be different
+%      20200325 implementing method for irregular grid, cleaning up
 ny = size(lon,1);
 nx = size(lon,2);
 np = numel(x);
 
-X = reshape(lon,[],1);
-Y = reshape(lat,[],1);
+% make sure entry points are column vectors
+x = x(:);
+y = y(:);
 
 % Get the grid spacing in x and y directions 
 % (trying both directions so could be meshgrid or ndgrid format)
-dx  = max(abs(lon(2,1)-lon(1,1)),abs(lon(1,2)-lon(1,1)));
-dy  = max(abs(lat(1,2)-lat(1,1)),abs(lat(2,1)-lat(1,1)));
+dx  = diff(lon(:,1));
+dy  = diff(lat(1,:));
 
-x = x(:);
-y = y(:);
+if max(dx) ~= min(dx) || max(dy) ~= min(dy)
+    % % IRREGULAR GRID (SLOWER)
+
+    % convert ndgrid to vector
+    lonV = reshape(lon(:,1),1,[]); 
+    % repeat the vector along the length of x
+    LON = repmat(lonV,np,1);
+    % find the closest lon to each point of x
+    [~,IX1] = min(abs(x - LON),[],2);
+
+    % convert ndgrid to vector
+    latV = lat(1,:); 
+    % repeat the vector along the length of y
+    LAT = repmat(latV,np,1);
+    % find the closest lon to each point of y
+    [~,IX2] = min(abs(y - LAT),[],2);
 
-IX1 = (x-X(1))./dx + 1;
-IX2 = (y-Y(1))./dy + 1;
+else
+    % % REGULAR GRID (FASTER)
+    dx = dx(1); dy = dy(1);
 
-IX1 = round(IX1);
-IX2 = round(IX2);
+    IX1 = (x-lon(1,1))/dx + 1;
+    IX2 = (y-lat(1,1))/dy + 1;
 
-IX1 = max([IX1,ones(np,1)],[],2);
-IX1 = min([IX1,ny*ones(np,1)],[],2);
-%
-IX2 = max([IX2,ones(np,1)],[],2);
-IX2 = min([IX2,nx*ones(np,1)],[],2);
+    IX1 = round(IX1);
+    IX2 = round(IX2);
+
+    IX1 = max([IX1,ones(np,1)],[],2);
+    IX1 = min([IX1,ny*ones(np,1)],[],2);
+
+    IX2 = max([IX2,ones(np,1)],[],2);
+    IX2 = min([IX2,nx*ones(np,1)],[],2);
+
+end
 
 IX = sub2ind([ny,nx],IX1,IX2);
-end
+
+end