sample_frequency=1;
%topo_file_name = 'SKOWHEGA.DEM';
%topo_file_name = 'CONCORD.DEM';
topo_file_name = 'VIDAL_CA-24000.dem';
%topo_file_name = 'MESQ.DEM';
%function [x_mat, y_mat, topo_mat, pixel_width] = show_24k_topo(topo_file_name, sample_frequency)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                     %%
%%  show_24k_topo.m                                    %%
%%                                                     %%
%%  This will plot a DEM format 1:24000 USGS topo map  %%
%%                                                     %%
%%  Arguments:                                         %%
%%    topo_file_name  : full name of the topo file     %%
%%    sample_frequency: take every n pixel             %%
%%                                                     %%
%%  Returned variables:                                %%
%%    x_mat           : x coordinate matrix (meters)   %%
%%    y_mat           : y coordinate matrix (meters)   %%
%%    topo_mat        : z coordinate matrix (meters)   %%
%%    pixel_width     : sample size of pixes (meters)  %%
%%                                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Declare variables  %%
%%%%%%%%%%%%%%%%%%%%%%%%%
pixel_spacing                  = 30;
pixel_width                    = sample_frequency * pixel_spacing;

%%%%%%%%%%%%%%%%%%%%%%%%
%%  Read in the data  %%
%%%%%%%%%%%%%%%%%%%%%%%%
[latgrat, longgrat, topo_mat, dem_info]   = usgs24kdem(topo_file_name, sample_frequency);

dem_info

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Announce the internal name of the DEM that we just read in  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(sprintf('Read DEM file: %s', getfield(dem_info, 'Quadranglename')))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Build some position matrices assuming the pixel spacing, %%
%%  flipping the y axis                                      %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[topo_mat_rows, topo_mat_cols]            = size(topo_mat);
x_vec                                     = sample_frequency * pixel_spacing * linspace(0, topo_mat_cols - 1, topo_mat_cols);
y_vec                                     = -sample_frequency * pixel_spacing * linspace(0, topo_mat_rows - 1, topo_mat_rows);
[x_mat, y_mat]                            = meshgrid(x_vec, y_vec);

%%%%%%%%%%%%%%%%%%%%%%%%
%%  Plot the results  %%
%%%%%%%%%%%%%%%%%%%%%%%%
%pcolor(x_mat, y_mat, topo_mat)
%surf(longgrat, -latgrat, topo_mat)
surf(x_mat, y_mat, topo_mat)
shading interp
daspect([1 1 1/5])
material shiny
camlight
xlabel('x [km]')
ylabel('y [km]')
zlabel('z [km]')
title(sprintf('%s, %d meter pixel spacing', getfield(dem_info, 'Quadranglename'), sample_frequency * pixel_spacing))
axis tight
view(0, 90)