diff --git a/README.md b/README.md
index f97e3301..ecd9013e 100644
--- a/README.md
+++ b/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
diff --git a/Tests/RunTests.m b/Tests/RunTests.m
index 8f9862b9..e98e2a07 100644
--- a/Tests/RunTests.m
+++ b/Tests/RunTests.m
@@ -15,6 +15,8 @@
    TestSanity 
 
    TestEleSizes
+   
+   TestInterp
 
    if exist('SRTM15+V2.1.nc')
    
diff --git a/Tests/TestInterp.m b/Tests/TestInterp.m
new file mode 100644
index 00000000..7110ba18
--- /dev/null
+++ b/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
diff --git a/datasets/CUDEMS/CUDEM_equal.nc b/datasets/CUDEMS/CUDEM_equal.nc
new file mode 100644
index 00000000..b45b0ab3
Binary files /dev/null and b/datasets/CUDEMS/CUDEM_equal.nc differ
diff --git a/datasets/CUDEMS/CUDEM_unequal.nc b/datasets/CUDEMS/CUDEM_unequal.nc
new file mode 100644
index 00000000..08d24ec3
Binary files /dev/null and b/datasets/CUDEMS/CUDEM_unequal.nc differ
diff --git a/datasets/CUDEMS/CUDEM_varying.nc b/datasets/CUDEMS/CUDEM_varying.nc
new file mode 100644
index 00000000..ef86c4a6
Binary files /dev/null and b/datasets/CUDEMS/CUDEM_varying.nc differ
diff --git a/utilities/FindLinearIdx.m b/utilities/FindLinearIdx.m
index b2bfe4da..e242c749 100644
--- a/utilities/FindLinearIdx.m
+++ b/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
\ No newline at end of file
+
+end