function d = getbathy(range, varargin)

% ***********************************************************************
% function d = getbathy(range, dflag)
%
% THIS FUNCTION EXTRACTS TOPOGRAPY DATA FROM THE DESIRED DATA SET
% NOTE: nctoolbox is required (http://code.google.com/p/nctoolbox/) &
% Setup_nctoolbox.m needs to be executed once during your Matlab session
%
% INPUT:
%   range = [lat1 lat2 lon1 lon2] = range array
%     lat1  = min latitude (-72 to 72)
%     lat2  = max lattitude (-72 to 72)
%     lon1  = west most longitude (degrees east or neg degrees west)
%     lon2  = east most longitude (degrees east or neg degrees west)
%
%   dflag = flag to choose data set (0=Smith & Sandwell, 1=etopo2, 2=arctic
%                                    3= etopo1, 4=S florida, 5=Monterey)
%           (if no flag the default will be 0=Smith & Sandwell)
%
% OUTPUT
%   d      = A structure of topgraphy data)
%   d.info = pertinent info, cell array
%   d.lat  = subset latitude (-90 to 90
%   d.lon  = subset longitude (degrees east)
%   d.topo = topography matrix, meters
%
% EXAMPLES (these are all equivalent)
%   NEPac = getbathy([35 65 -160 -120], 0);
%   NEPac = getbathy([35 65 200 240], 0);
%   NEPac = getbathy([35 65 200 240]);
%
% Check with: 
% contourf(NEPac.lon,NEPac.lat,NEPac.topo,'linecolor','none');
% ***********************************************************************

% ***********************************************************************
% ***********************************************************************
%   DEFINE DATA PATHES (ADD OR CHANGE AS YOU SEE FIT)
% ************************************************************************
topo_path_0 = 'http://ferret.pmel.noaa.gov/thredds/dodsC/data/PMEL/smith_sandwell_topo_v8_2.nc';
topo_path_1 = 'http://thredds.socib.es/thredds/dodsC/bathymetries/etopo2.nc';
topo_path_2 = ['C:\Documents and Settings\jplant\My Documents\',...
    'MATLAB\PWP MODEL\Arctic7564\IBACO\IBCAO_ver2_23_GEO_netCDF_2min.grd'];
topo_path_3 = 'http://ourocean.jpl.nasa.gov:8080/thredds/dodsC/hycom/etopo1.nc';
%topo_path_4 = 'C:\Users\jplant\Documents\MATLAB\PWP MODEL\Maps\data\S_florida.nc';
topo_path_4 = 'C:\Users\jplant\Documents\MATLAB\PWP MODEL\Maps\data\Tampa.nc';
topo_path_5 = 'C:\Users\jplant\Documents\MATLAB\Tuna\data\Monterey.nc';


%topo_path_0 = 'http://las.pfeg.noaa.gov/thredds/dodsC/bathy/smith_sandwell_topo.nc';
%topo_path_1 = 'http://las.pfeg.noaa.gov/thredds/dodsC/bathy/etopo2.nc';

% ************************************************************************
%   CHECK FOR DATA SET FLAG ( DEFAULT = 0 = Smith & Sandwell topography)
% ************************************************************************
if size(varargin,2) == 0
    dflag = 0;
elseif size(varargin,2) ==1
    dflag = varargin{1};
else
    disp('Data flag should be a scalar')
    return
end

% ************************************************************************
% CHECK LONGITUDE: DEGREES EAST? ZERO MERIDIAN CROSSED -180/+180 CROSSED
% ************************************************************************
% dflag = 4; % for testing
% bnds  = [24.7 26 -81.6 -80.0]; % for testing

bnds = [range(1) range(2) range(3) range(4)];

if dflag == 0 % 0 to 360 E
    T1 = bnds(3:4) < 0;     % Look for negative values in longitude
    bnds(3:4) = bnds(3:4) + T1*360; % Convert to 0 to 360 degrees east
else % -180 to + 180
    T1 = bnds(3:4) > 180;     % Look for values in longitude > +180
    bnds(3:4) = bnds(3:4) - T1*360; % Convert to -180 to +180 degrees east
end


lon_flag = 0;           % Default no meridian or  +180 to -180 crossed
if bnds(3) > bnds(4)    % 0 meridian or +180 to -180 crossed !!!!
    lon_flag = 1;
end

% ************************************************************************
%           CHOOSE DATA SET BASED ON FLAG VALUE ( DEFAULT = 0)
% ************************************************************************
if dflag == 0 % Smith & Sandwell
    BATHYpath = topo_path_0;
    ds  = ncdataset(BATHYpath); % create data ID
    lon = double(ds.data('longitude')); % get longitude (degrees east) 0 to 360
    lat = double(ds.data('latitude')); % get latitude 
    dname='ROSE';
    disp('Connecting to Smith & Sandwell topography');
    disp(BATHYpath);
elseif dflag == 1 % etopo2
    BATHYpath = topo_path_1;
    ds  = ncdataset(BATHYpath); % create data ID
    lon = double(ds.data('lon')); % get longitude (degrees east) -180 to +180
    lat = double(ds.data('lat')); % get latitude 
    dname ='topo';
    disp('Connecting to etopo2 topography');
    disp(BATHYpath);  
elseif dflag == 2 % Arctic topo (IBCAO V2.23 2min)
    BATHYpath = topo_path_2;
    ds      = ncdataset(BATHYpath); % create data ID    
    x_range = ds.data('x_range');
    y_range = ds.data('y_range');
    dim     = ds.data('dimension');
    
    lon     = linspace(x_range(1), x_range(2),dim(1)); %-180 to +180
    lat     = linspace(y_range(2), y_range(1),dim(2))'; % 90 to 54
    dname   = 'z';
    disp('Connecting to Arctic topography (IBCAO V2.23 2min)');
    disp(BATHYpath); 
elseif dflag == 3 % Setopo1
    BATHYpath = topo_path_3;
    ds  = ncdataset(BATHYpath); % create data ID
    lon = double(ds.data('lon')); % get longitude (degrees east) -180 to +180
    lat = double(ds.data('lat')); % get latitude 
    dname ='zz';
    disp('Connecting to etopo1 topography');
    disp(BATHYpath);  
elseif dflag == 4 % S. Florida
% The affine transform consists of six coefficients returned by 
% GDALDataset::GetGeoTransform() which map pixel/line coordinates into
% georeferenced space using the following relationship:

%    Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
%    Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

% In case of north up images, the GT(2) and GT(4) coefficients are zero, and
% the GT(1) is pixel width, and GT(5) is pixel height. The (GT(0),GT(3))
% position is the top left corner of the top left pixel of the raster.
% Note that the pixel/line coordinates in the above are from (0.0,0.0)
% at the top left corner of the top left pixel to (width_in_pixels,height_in_pixels)
% at the bottom right corner of the bottom right pixel. The pixel/line
% location of the center of the top left pixel would therefore be (0.5,0.5).

    BATHYpath = topo_path_4;
    ds    = ncdataset(BATHYpath); % create data ID
    dname = 'Band1'; % topo grid
    yx  = ds.size(dname); % row , col or y x
    tmp = ds.attributes('GDAL_Geographics');
    geotrans = tmp{find(strcmp('GeoTransform',tmp(:,1))),2};
    geotrans = str2num(geotrans);
    lon = double(0:yx(2)-1) * geotrans(2) + geotrans(1) + geotrans(2)/2; %pixel center
    lat = double(0:yx(1)-1) * geotrans(6) + geotrans(4) + geotrans(6)/2; %pixel center
    disp('Connecting to S. Florida subset topography');
    disp(BATHYpath);  
elseif dflag == 5 % Monterey
    BATHYpath = topo_path_5;
    ds    = ncdataset(BATHYpath); % create data ID
    dname = 'Band1'; % topo grid
    yx  = ds.size(dname); % row , col or y x
    tmp = ds.attributes('GDAL_Geographics');
    geotrans = tmp{find(strcmp('GeoTransform',tmp(:,1))),2};
    geotrans = str2num(geotrans);
    lon = double(0:yx(2)-1) * geotrans(2) + geotrans(1) + geotrans(2)/2; %pixel center
    lat = double(0:yx(1)-1) * geotrans(6) + geotrans(4) + geotrans(6)/2; %pixel center
    disp('Connecting to Monterey Bay subset topography');
    disp(BATHYpath);  
else
    disp('Data flag not recognized: only 0 or 1 at this point.')
    return
end
info{1,1} =BATHYpath;

% ************************************************************************
%           GET INDICES FOR SUBSET & EXTRACT DATA
% ************************************************************************
lon_ind1 = find(lon < bnds(3),1,'last');
lon_ind2 = find(lon > bnds(4),1,'first');

if lat(1) < lat(end) %low to high
    lat_ind1 = find(lat < bnds(1),1,'last');
    lat_ind2 = find(lat > bnds(2),1,'first');
else %high to low (North Up)
    lat_ind1 = find(lat >= bnds(2),1,'last');
    lat_ind2 = find(lat <= bnds(1),1,'first');
end

if dflag ~= 2
    d.lat  = lat(lat_ind1:lat_ind2); % latitude subset
    if lon_flag ==0 % NO 0 MERIDIAN CROSSING
        d.lon  = lon(lon_ind1:lon_ind2); % longitude subset
        d.topo = double(ds.data(dname,[lat_ind1 lon_ind1],[lat_ind2 lon_ind2]));
    else            % 0 MERIDIAN CROSSED!!!!!
        d.lon = [lon(lon_ind1:length(lon));lon(1:lon_ind2)];
        d1 = double(ds.data(dname,[lat_ind1 lon_ind1],[lat_ind2 length(lon)])); % max lon to 360
        d2 = double(ds.data(dname,[lat_ind1 1],[lat_ind2 lon_ind2])); % 0 ti min lon
        d.topo = cat(2,d1,d2);
        disp('Data range crosses zero meridian.')
    end
% Arctic netcdf a little different - need to parse out the  1 d array with
% a loop to get it into a matrix   
else 
    d.lat  = lat(lat_ind1:lat_ind2); % latitude subset
  
   if lon_flag ==0 % NO 0 MERIDIAN CROSSING 
       d.lon  = lon(lon_ind1:lon_ind2); % longitude subset
       temp = ones(length(d.lat),length(d.lon)) * NaN; % predimmension
       ct=1;
       for i = 1: (lat_ind2-lat_ind1)+1
           ind1 = (lat_ind1+i-1) * dim(1) +1;
           ind2 = (lat_ind1 + i) * dim(1);
           tmp  = double(ds.data(dname,ind1,ind2)); % extract row
           temp(i,:) = tmp(lon_ind1:lon_ind2);
       end
       d.topo =temp;
   else
      % STILL NEED TO ADD THE +180 to -180  EXTRACTION CODE
       
       
   end
end





