% paleoflood.m
% Surveyed data from Queen Creek, 10/15/00 Geomorphology
%JRA, 10/16/00
load flot.txt
load bottom.txt
load slope.txt
load xsec.txt

dn = 2.5;
de = dn;

figure(1)
clf
%Flotsam
plot(flot(:,2),flot(:,3),'k*');
hold on
for i = 1:length(flot(:,1))
   s = sprintf('%.f', flot(i,1));
   text(flot(i,2)+de,flot(i,3)+dn,s);
   drawnow; 
end

%bottom
plot(bottom(:,2),bottom(:,3),'k^');
for i = 1:length(bottom(:,1))
   s = sprintf('%.f', bottom(i,1));
   text(bottom(i,2)+de,bottom(i,3)+dn,s);
   drawnow; 
end

%slope
plot(slope(:,2),slope(:,3),'kv');
for i = 1:length(slope(:,1))
   s = sprintf('%.f', slope(i,1));
   text(slope(i,2)+de,slope(i,3)+dn,s);
   drawnow; 
end

%xsec
plot(xsec(:,2),xsec(:,3),'k>');
for i = 1:length(xsec(:,1))
   s = sprintf('%.f', xsec(i,1));
   text(xsec(i,2)+de,xsec(i,3)+dn,s);
   drawnow; 
end

axis([800 990 860 1060])
axis equal
ylabel('Northing (m)')
xlabel('Easting (m)')
title('Queen Creek Paleo flood survey data')
legend('Flotsam','Channel bottom','Channel slope','Channel x-section',4)

print -depsc points.eps



%data for the gridding.  We want to have the same extents for the grid so that it can be subtracted.
gdata =  [bottom; xsec; slope; flot];
ge=gdata(:,2);
gn=gdata(:,3);


%produce contour map of channel
figure(2)
clf
%The loaded file is in the array  NETR
data =  [bottom; xsec; slope];
%refer to columns as:
%data(:,1) is point number
%data(:,2) is Easting
%data(:,3) is Northing
%data(:,4) is Elevation

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

%AND HERE
e=data(:,2);
n=data(:,3);
h=data(:,4);
estep = abs((min(ge)-max(ge))/nx);
nstep = abs((min(gn)-max(gn))/ny);
ee = [min(ge):estep:max(ge)];
nn = [min(gn):nstep:max(gn)];
[xi,yi,zi] = griddata(e,n,h,ee,nn');
chgrid=zi;
v = floor(min(h)):ci:max(h);

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

%label contours
clabel(c,h);
colormap([0,0,0])
hold on

%show data:
plot(e,n,'<')
axis([800 990 860 1060])
xlabel('easting (m)')    
ylabel('northing (m)')   
title('Elevation and observation points--CHANNEL')
axis equal

print -depsc channel.eps


%produce contour map of flotsam
figure(3)
clf
%The loaded file is in the array  NETR
data = flot;
%refer to columns as:
%data(:,1) is point number
%data(:,2) is Easting
%data(:,3) is Northing
%data(:,4) is Elevation

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

%AND HERE
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');
flotgrid=zi;

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

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

%label contours
clabel(c,h);
colormap([0,0,0])
hold on

%show data:
plot(e,n,'*')
axis([800 990 860 1060])
xlabel('easting (m)')    
ylabel('northing (m)')   
title('Elevation and observation points--FLOTSAM')
axis equal

print -depsc flot.eps

%Subtract the grids
figure(4)
clf
plot(900,900,'k.')
hold on
depth=flotgrid-chgrid;
imagesc(ee,nn',depth)
colormap(gray)
colorbar
plot(e,n,'*')
axis([800 990 860 1060])
xlabel('easting (m)')    
ylabel('northing (m)')   
title('Depth: FLOTSAM-CHANNEL (meters)')
axis equal

print -depsc depth.eps