%script to produce several different plots of the Van Matire topographic data  in an illustration of contouring.
%JRA August 26, 2001

load -ascii vm.txt

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

%CHANGE THINGS BETWEEN HERE
%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 = 100;
ny = nx;
%contour interval:
ci = 2; 

%Check this assignment for the columns in the loaded datafile.
n=data(:,2);
e=data(:,3);
hh=data(:,4);
%AND HERE



%Here is the map of points:
figure(1)
clf
plot(e, n, 'k+')
xlabel('easting')    
ylabel('northing')   
title('Observation points')
axis equal
print -depsc vmtopopts.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')    
ylabel('northing')   
title('Grid and observation points')
axis equal

print -depsc vmptsandgrid.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

print -depsc vmcontourslin.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
print -depsc vmtriangles.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

print -depsc vmcontourscubic.eps


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

xlabel('easting')    
ylabel('northing')   
title('Cubic interpolation and observation points')
axis equal
print -depsc vmimage.eps  


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