close all
fprintf(1,'Reading Arcgrid ...\n');
%[Z,R] = arcgridread('montbays15asc.asc');
Zsave = Z;
 
%remove the rows and columns with all NaN
fprintf(1, 'Removing NaN rows and columns ...\n');
a = isfinite(Zsave);

%find columns with any real finite numbers
c = any(a,2);
%find first column with a finite number
cFirst = min(find(c==1));
%find last column with a finite number
cLast = max(find(c==1)); 

%find rows with any real finite numbers
r = any(a,1);
%find first row with a finite number
rFirst = min(find(r==1));
%find last column with a finite number
rLast = max(find(r==1)); 

Z2 = Zsave(cFirst:cLast, rFirst:rLast);

%Re-map R array
fprintf(1, 'Generating new Ref array ...\n');
[x11 y11] = pix2map(R,cFirst,rFirst);
R2 = makerefmat(x11, y11,R(1:2,1),R(1:2,2));

fprintf(1, 'Getting Haxby color map ...\n');
cm = haxby(256);

fprintf(1, 'Showing map ...\n');
mapshow(Z2,R2,'DisplayType','surface');
xlabel('x (easting in meters)'); 
ylabel('y (northing in meters)')
colormap(cm);

% Format the map
fprintf(1, 'Formating the map ...\n');
axis normal; 
view(2); 
axis equal; 
grid on
zlabel('elevation in meters')

%Decimate the map
fprintf(1, 'Decimating the map ...\n');
Z2decimated = resizem(Z,0.1);
figure
mapshow(Z2decimated,R2,'DisplayType','surface');
xlabel('x (easting in meters)'); 
ylabel('y (northing in meters)')
colormap(cm);

% View the terrain in 3-D
axis normal; 
view(2); 
axis equal; 
grid on
zlabel('elevation in meters')
