function hdf_struct = load_hdf_coords( hdf_filename )

% This function was adapted from 'load_hdf_file.m' by Dave Morrish of
% University of Western Ontario, which was distributed as part of the 
% NIST Database Software. The portion for importing time series was 
% removed, because time series are not needed for assembling the
% tributary area matrix.

if ~exist('hdf_filename') || isempty(hdf_filename)
    % Load pathname for HDF files used in last analysis, if available, otherwise create empty string:
    if exist('tmp_HDF_path.mat','file'), load('tmp_HDF_path'); end
    if ~exist('HDF_path','var') || ~isstr(HDF_path), HDF_path = ''; end
    % Prompt use to select pressure database file:
    [HDF_filename HDF_path] = uigetfile({'*.hdf','HDF files (*.hdf)'},'Select an HDF pressure database file in the directory of interest',HDF_path);
    if any(HDF_filename~=0)
        hdf_filename = fullfile(HDF_path, HDF_filename);
        save('tmp_HDF_path','HDF_path');
    else
        err = errordlg('HDF pressure database file selection cancelled', 'Analysis cancelled');
        uiwait(err); hdf_struct.error = 1; return;
    end
end

hdf_struct.File_Name=[];
sd_id=hdfsd('start',hdf_filename,'read'); % open hdf file
if sd_id <=0
    error(['Could not open >',hdf_filename,'< for reading: check file name, paths, and permisions.']);
end
[ndatasets,natt,sata]=hdfsd('fileinfo',sd_id); % get file info
for sds_idx=0:ndatasets-1
    sds_id=hdfsd('select',sd_id,sds_idx);
    [dsname,dsndim,dsdims,dstype,dsatts,stat]=...
        hdfsd('getinfo',sds_id);
    ds_start=zeros(1,dsndim);
    ds_stride=[];
    ds_edges=dsdims;
    if ~strcmp(dsname,'Time_Series') % do not output time series data
        [ds_data,status]=hdfsd('readdata',sds_id,ds_start,ds_stride,ds_edges);
        %
        % retrive attributes
        %
        for i=1:dsatts
            [aname,adt,count,status]=hdfsd('attrinfo',sds_id,i-1);
            [attrib,status]=hdfsd('readattr',sds_id,i-1);
            aname=strrep(aname,' ','_');
            if status==0
                hdf_struct=setfield(hdf_struct,['A_',dsname,'_',aname],attrib);
            end
        end
            hdf_struct=setfield(hdf_struct,dsname,ds_data);
    end
end
hdfml('closeall');

% Correct known errors in tap coordinates:
if strncmp(hdf_struct.File_Name,'ADW100o100D024a',15)
    % 3D coordinates for tap 1703 are missing:
    if isempty(find(hdf_struct.Tap_Coordinates_3D(:,1)==1703));
        hdf_struct.Tap_Coordinates_3D = [hdf_struct.Tap_Coordinates_3D; 1703 5 11.64 74.5 24.97];
    end
elseif strncmp(hdf_struct.File_Name, 'ADW600o100D024a',15)
    % y-coordinates for taps 516 and 616 are 130.75 but should be 124.5:
    hdf_struct.Tap_Coordinates_3D(find(hdf_struct.Tap_Coordinates_3D(:,1)==516 | hdf_struct.Tap_Coordinates_3D(:,1)==616),4) = 124.5;
    % y-coordinates for taps 916 and 1016 are 105.75 but should be 112:
    hdf_struct.Tap_Coordinates_3D(find(hdf_struct.Tap_Coordinates_3D(:,1)==916 | hdf_struct.Tap_Coordinates_3D(:,1)==1016),4) = 112;
    % y-coordinates for taps 1116 and 1216 are 93.25 but should be 105.75:
    hdf_struct.Tap_Coordinates_3D(find(hdf_struct.Tap_Coordinates_3D(:,1)==1116 | hdf_struct.Tap_Coordinates_3D(:,1)==1216),4) = 105.75;
    % y-coordinates for taps 1316 and 1401 are 80.75 but should be 99.5:
    hdf_struct.Tap_Coordinates_3D(find(hdf_struct.Tap_Coordinates_3D(:,1)==1316 | hdf_struct.Tap_Coordinates_3D(:,1)==1401),4) = 99.5;
    % y-coordinates for taps 1406 and 1407 are 112 but should be 118.25:
    hdf_struct.Tap_Coordinates_3D(find(hdf_struct.Tap_Coordinates_3D(:,1)==1406 | hdf_struct.Tap_Coordinates_3D(:,1)==1407),4) = 118.25;
    % y-coordinates for taps 1410 and 1411 are 99.5 but should be 112:
    hdf_struct.Tap_Coordinates_3D(find(hdf_struct.Tap_Coordinates_3D(:,1)==1410 | hdf_struct.Tap_Coordinates_3D(:,1)==1411),4) = 112;
    % y-coordinates for taps 1414 and 1415 are 87 but should be 105.75:
    hdf_struct.Tap_Coordinates_3D(find(hdf_struct.Tap_Coordinates_3D(:,1)==1414 | hdf_struct.Tap_Coordinates_3D(:,1)==1415),4) = 105.75;

end