%script to import a text file 
%and make a contour map.  It will not as of yet do breaklines.

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

%The loaded file is in the array  NETR
data = ermatxyz;
%refer to columns as:
%data(:,1) is block number
%data(:,2) is 116
%data(:,3) is point number
%data(:,4) is Easting
%data(:,5) is Northing
%data(:,6) is Elevation
%data(:,7) is code

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

%AND HERE





%Here is the map of points:
figure(1)
clf
plot(data(:,2), data(:,3), 'k+')

%Here is contour map:
figure(2)
clf
e=data(:,2);
n=data(:,3);
h=data(:,4);
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,h,ee,nn');

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


[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 (m)')    
ylabel('northing (m)')   
title('Elevation and observation points')
axis equal





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

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

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

xlabel('easting (m)')    
ylabel('northing (m)')   
title('Elevation and observation points')
axis equal
axis([875 1150 825 1125])


figure(3)
%now the trench outlines
load trenchN.txt
load trenchS.txt
plot(trenchS(:,1),trenchS(:,2), 'r-')
plot(trenchN(:,1),trenchN(:,2), 'r-')

%now the wall at the base of the scarp
load wall1.txt
load wall2.txt
load wall3.txt
plot(wall1(:,1),wall1(:,2), 'k-')
plot(wall2(:,1),wall2(:,2), 'k-')
plot(wall3(:,1),wall3(:,2), 'k-')



figure(4)
clf
surfl(xi,yi,zi)
shading interp
view(54,44)
colormap pink
hold on
plot3(trenchS(:,1),trenchS(:,2),trenchS(:,3)+.2)
plot3(trenchN(:,1),trenchN(:,2),trenchN(:,3))
%plot3(data(:,2),data(:,3),data(:,4), 'k.')
xlabel('Easting (m)');ylabel('Northing');zlabel('elevation (m)')


%export this one
figure(3)
print -depsc contoursfrommatlab.eps


%figure(5)
%plot(1000,1000)
%imagesc(ee,nn,zi')
%axis equal