%script to produce several different plots of the DAVIS data from the Middleton book in an illustration of contouring.
%this is mostly black and white.

load -ascii DMAP.DAT

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

%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 = 20;
ny = nx;
%contour interval:
ci = 20; 

%AND HERE

n=data(:,1);
e=data(:,2);
hh=data(:,3);

%Here is the map of points:
figure(1)
clf
plot(e, n, 'k+')
xlabel('easting')    
ylabel('northing')   
title('Observation points')
axis equal
axis([0 6.5 0 6.5])
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);

hold on

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

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

xlabel('easting')    
ylabel('northing')   
title('Grid and observation points')
axis equal
axis([0 6.5 0 6.5])
colormap([0 0 0])
print -depsc ptsandgrid.eps


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

%label contours
clabel(c,h);
%use this one if you want to label the contours manually:
%clabel(c,h, 'manual');
colormap([0  0 0])
hold on

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

xlabel('easting')    
ylabel('northing')   
title('Linear interpolation and observation points')
axis equal
axis([0 6.5 0 6.5])
print -depsc contourslin.eps

%triangles for the linear interpolation
figure(4)
clf
tri = delaunay(e,n);
trimesh(tri,e,n,hh)
colormap([0 0 0])
view(2)
xlabel('easting')    
ylabel('northing')   
title('Triangles for linear interpolation')
axis equal
axis([0 6.5 0 6.5])
print -depsc triangles.eps

%nice contour map
figure(5)
clf
[xi,yi,zi] = griddata(e,n,hh,ee,nn', 'cubic');
[c,h] = contour(xi,yi,zi,v);

%label contours
clabel(c,h);
%use this one if you want to label the contours manually:
%clabel(c,h, 'manual');
colormap([0 0 0])
hold on

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

xlabel('easting')    
ylabel('northing')   
title('Cubic interpolation and observation points')
axis equal
axis([0 6.5 0 6.5])
print -depsc contourscubic.eps


%imagesc
figure(6)
clf
plot(0,0,'k+')
hold on
imagesc(ee',nn',zi)
colormap(gray)
%show data:
plot(e,n,'k.')

xlabel('easting')    
ylabel('northing')   
title('Cubic interpolation and observation points')
axis equal
axis([0 6.5 0 6.5])
print -depsc image.eps  


%surface
figure(7)
clf
surfl(ee,nn,zi)
shading interp
colormap(gray)
hold on
plot3(e,n,hh,'b.')
view(116,48)
xlabel('easting')    
ylabel('northing')
zlabel('elevation')   
title('Cubic interpolation and observation points')
%print -depsc surface.eps 

%show the points with the data values
figure(8)
clf
plot(e,n, 'k.')
hold on
for i=1:length(e)
     s = sprintf('%d\n',hh(i));
     text(e(i),n(i),s)
end
xlabel('easting')    
ylabel('northing')
title('Observation points and data values')
axis([0 6.5 0 6.5])
print -depsc mapwithpointslabeled.eps
