Code covered by the BSD License

# Surface2VTK

### James Ramm (view profile)

Interpolates scattered data and saves the result to a vtk file.

fillElevation(Elev,ncellX,ncellY,nlayer,nR)
```function [Elev]= fillElevation(Elev,ncellX,ncellY,nlayer,nR)

%==========================================================================
% All values of nR must have 4 surrounding elevation points with real
% values. This function runs through relevant points and replaces any
% occuring NaN values with a mean of the surrounding points.
%==========================================================================

%% Check grid corner points
for l = 1:nlayer+1
if isnan(Elev(1,1,l))
Elev(1,1,l)=mean([Elev(1,2,l),Elev(2,1,l),Elev(2,2,l)]);

elseif isnan(Elev(1,ncellY+1,l))
Elev(1,ncellY+1,l) = mean([Elev(1,ncellY,l),Elev(2,ncellY+1,l),Elev(2,ncellY,l)]);

elseif isnan(Elev(ncellX+1,1,l))
Elev(ncellX+1,1,l) = mean([Elev(ncellX,1,l),Elev(ncellX+1,2,l),Elev(ncellX,2,l)]);

elseif isnan(Elev(ncellX+1,ncellY+1,l))
Elev(ncellX+1,ncellY+1,l)=mean([Elev(ncellX+1,ncellY,l),Elev(ncellX,ncellY+1,l),Elev(ncellX,ncellY,l)]);
end
end
%% Check inner points
fprintf('assigning missing elevation points to cell vertices...       \n\n');
for j = 2:ncellY
for i = 2:ncellX
for l = 1:nlayer+1
if isfinite(nR(i,j,1))
if isnan(Elev(i,j,l))
Elev(i,j,l) = mean([Elev(i-1,j,l),Elev(i+1,j,l),Elev(i,j-1,l),Elev(i,j+1,l)]);
end
end
end
end
end

%% Check top and bottom edges (1st and last row)
for j = 2:ncellY
for l=1:nlayer+1
if isnan(Elev(1,j,l))
Elev(1,j,l)=mean([Elev(1,j+1,l),Elev(1,j-1,l),Elev(2,j,l)]);
end
if isnan(Elev(ncellX+1,j))
Elev(ncellX+1,j,l) = mean([Elev(ncellX,j,l),Elev(ncellX+1,j+1,l),Elev(ncellX+1,j-1,l)]);
end
end
end

%% Check sides (first and last column)
fprintf ('adjusting E-W boundary vertices...      \n\n');
for i=2:ncellX
for l = 1:nlayer+1
if isnan(Elev(i,1,l))
Elev(i,1,l)=mean([Elev(i,2,l),Elev(i-1,1,l),Elev(i+1,1,l)]);
end
if isnan(Elev(i,ncellY+1,l))
Elev(i,ncellY+1,l)=mean([Elev(i,ncellY),Elev(i-1,ncellY+1),Elev(i+1,ncellY+1,l)]);
end
end
end

%% Check for elevation overlaps (may occur as a result of interpolation)
corELEV=0;
for l=2:nlayer+1
fprintf('correcting elevation overlaps in layer %i.....\n',l)
for j = 1:ncellY+1
for i = 1:ncellX+1
if isfinite(Elev(i,j,l));
if Elev(i,j,l)>Elev(i,j,l-1)
Elev(i,j,l)=Elev(i,j,l-1);
corELEV = corELEV+1;
end
end
end
end
fprintf('%i Elevation points have been corrected in layer %i\n',corELEV,l);
end
%end of function
end

```