function [data,lonax,latax] = readetopo1(lon1, lat1, lon2, lat2, p)
% Hi! READETOPO1 lets you read data from the Etopo1 global relief database.
%
% Usage:
% [DATA, LONAX, LATAX] = READETOPO1(LON1, LAT1, LON2, LAT2, DRAWMAP)
%
% LON1 and LAT1 are the upper left coordinates of the area of interest, and
% LON2 and LAT2 are the lower right coordinates of the area of interest.
% Presently the function cannot handle crossing the limits of the map. Set
% the switch DRAWMAP to a nonzero value to get a plot of the topography and
% bathymetry.
%
% Example:
% [DATA, LONAX, LATAX] = READETOPO1(3, 54, 6, 51, 1) will give data, axes
% and a plot of the area around the Netherlands.
% Use fractions if you want, but it is also possible to enter degrees and
% minutes of arc as a vector; 55.5 is the same as [55 30], for example.
%
% Copywrong this script Jeroen Janmaat, 2010.
% First version 2010.09.23.
%
% The Etopo1 database you need is the BINARY and GRID/NODE-REGISTERED
% version that can be downloaded from the NOAA website.
% "Amante, C. and B. W. Eakins, ETOPO1 1 Arc-Minute Global Relief Model:
% Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS
% NGDC-24, 19 pp, March 2009."

filename = 'etopo1_bed_g_i2.bin';
griddims = [21601 10801];
samplesize = 2;

lon1 = convlon1(lon1);
lon2 = convlon2(lon2);
lat1 = convlat1(lat1);
lat2 = convlat2(lat2);

if(lon2(1) + lon2(2)/60) < (lon1(1) + lon1(2)/60)
    disp('boe!')
end

% Open datafile.
fid = fopen(filename,'r');

% Skip lines of data.
fseek(fid,griddims(1)*samplesize*(lat1(1)*60 + lat1(2)),'bof');

% Skip to leftmost longitude.
fseek(fid,samplesize*(lon1(1)*60+lon1(2)),'cof');

% Define data matrix
data = zeros(((lat2(1)*60 + lat2(2)) - (lat1(1)*60 + lat1(2))) + 1, ...
    (lon2(1)*60 + lon2(2)) - (lon1(1)*60 + lon1(2)) + 1);

% Now, read data
for ii = 1:((lat2(1)*60 + lat2(2)) - (lat1(1)*60 + lat1(2))) + 1
    data(ii,:) = fread(fid,(lon2(1)*60 + lon2(2)) - (lon1(1)*60 + lon1(2)) + 1,'int16');
    fseek(fid,samplesize*(griddims(1)-((lon2(1)*60 + lon2(2)) - (lon1(1)*60 + lon1(2))) - 1),'cof');
end

% Close datafile.
fclose(fid);

rlon1 = single((lon1(1)+lon1(2)/60) - 180);
rlat1 = single(90 - (lat1(1)+lat1(2)/60));
rlon2 = single((lon2(1)+lon2(2)/60) - 180);
rlat2 = single(90 - (lat2(1)+lat2(2)/60));
lonax = rlon1:1/60:rlon2;
latax = rlat1:-1/60:rlat2;

% Plot if switch 'p' is nonzero
if p
    figure
    imagesc(lonax,latax,data)
    axis xy
    xlabel('Longitude [degrees]')
    ylabel('Latitude [degrees]')
    set(gca,'DataAspectRatio',[1 cosd(mean(latax)) 1])
end

end









function lon1 = convlon1(lon1)
% Straighten out leftmost longitude. Convert to pseudo-longitude.
if length(lon1) == 1
    lon1(1) = lon1(1) + 180;
    lon1(2) = floor((lon1 - single(floor(lon1)))*60);
    lon1(1) = floor(lon1(1));
else
    if lon1(1) < 0
        lon1(1) = lon1(1) + 179;
        lon1(2) = 60 - lon1(2);
    else
        lon1(1) = lon1(1) + 180;
    end
end
end

function lon2 = convlon2(lon2)
% Straighten out rightmost longitude. Convert to pseudo-longitude.
if length(lon2) == 1
    lon2(1) = lon2(1) + 180;
    lon2(2) = ceil((lon2 - single(floor(lon2)))*60);
    lon2(1) = floor(lon2(1));
    if lon2(2) == 60
        lon2(2) = 0;
        lon2(1) = lon2(1) + 1;
    end
else
    if lon2(1) < 0
        lon2(1) = lon2(1) + 179;
        lon2(2) = 60 - lon2(2);
    else
        lon2(1) = lon2(1) + 180;
    end
end
end

function lat1 = convlat1(lat1)
% Straighten out upper latitude. Convert to pseudo-latitude.
if length(lat1) == 1
    lat1(1) = 90 - lat1(1);
    lat1(2) = floor((lat1 - single(floor(lat1)))*60);
    lat1(1) = floor(lat1(1));
else
    if lat1(1) < 0
        lat1(1) = 90 - lat1(1);
    else
        lat1(1) = 89 - lat1(1);
        lat1(2) = 60 - lat1(2);
    end
end
end

function lat2 = convlat2(lat2)
% Straighten out lower latitude. Convert to pseudo-latitude.
if length(lat2) == 1
    lat2(1) = 90 - lat2(1);
    lat2(2) = ceil((lat2 - single(floor(lat2)))*60);
    lat2(1) = floor(lat2(1));
    if lat2(2) == 60
        lat2(2) = 0;
        lat2(1) = lat2(1) + 1;
    end
else
    if lat2(1) < 0
        lat2(1) = 90 - lat2(1);
    else
        lat2(1) = 89 - lat2(1);
        lat2(2) = 60 - lat2(2);
    end
end
end
