%script to take scattered topographic data
%and make a contour map.  It will not as of yet do breaklines.
%JRA, August 23, 2001

%CHANGE things between here:
%file to load is NETR.dat.txt or whatever
load -ascii filtered.txt

%The loaded file is in the array 
data = filtered;

%set parameters:
%number of x and y nodes for grid
%make these smaller if you get an error about excceding the max array size
nx = 200;
ny = nx;
%contour interval:
ci = 2; %feet

%AND HERE
%l=length(data);
%e=data(1:l,2)-min(data(:,2));
%n=data(1:l,3)-min(data(:,3));
%hh=data(1:l,4)-min(data(:,4));
n=data(:,2);
e=data(:,3);
hh=data(:,4);

%Here is the map of points:
figure(1)
clf
plot(e, n, 'k+')
xlabel('easting (ft)')    
ylabel('northing (ft)')   
title('Observation points')
axis equal
print -depsc topopts.eps

%Here is contour map:
figure(2)
clf
estep = abs((min(e)-max(e))/nx);
nstep = abs((min(n)-max(n))/ny);
ee = [min(e):estep:max(e)];
nn = [min(n):nstep:max(n)];
[xi,yi,zi] = griddata(e,n,hh,ee,nn');

v = floor(min(hh)):ci:max(hh);


[c,h] = contour(xi,yi,zi,v);

%label contours
clabel(c,h);
colormap(jet)
hold on

%show data:
plot(e,n,'o')

%show grid:
plot(xi,yi,'+')

xlabel('easting (ft)')    
ylabel('northing (ft)')   
title('Elevation and observation points')
axis equal

print -depsc ptsandgrid.eps




%nice contour map
figure(3)
clf
[c,h] = contour(xi,yi,zi,v);

%label contours
clabel(c,h, 'manual');
colormap(jet)
hold on

%show data:
%plot(e,n,'.')

xlabel('easting (ft)')    
ylabel('northing (ft)')   
title('Elevation and observation points, 2 foot contour interval')
axis equal

print -depsc contours.eps


%imagesc
figure(4)
clf
plot(0,0,'k+')
hold on
%imagesc(ee',nn',zi)
imagesc(zi)
axis equal
print -depsc image.eps  


%surface
figure(5)
clf
surfl(ee,nn,zi)
shading interp
colormap(gray)
hold on
%plot3(e,n,hh,'b.')