Code covered by the BSD License  

Highlights from
Surface2VTK

Surface2VTK

by

 

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
fprintf('adjusting corner vertices...     \n');
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)
fprintf('adjusting N-S boundary vertices...      \n\n'); 
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


    
    

Contact us