   function surGrid(sbeFile, column)
% surGrid - Create a contour of data along the CA coast from Pt Sur underway data
%
% Use as:  surGrid(dataFile, column)
%
% Inputs: dataFile = name of the underway data file
%         column      = Column of data in the seaBird data file to grid
%         dataFile = An underway data file containing two columns of data. The first
%                    column is time (in Excel's format - serial days since 1900). The second
%                    column is the data points.
%
% Example: surGrid('example_tsg_61200.txt', 6)
%          surGrid('Example_seabird.prc', 'Example_pco2.prn')
%%  Key ocean surface data columns:
%%	30 = SST
%%	31 = Fluor volts 
%%	33 = Salinity
%
% Note the graph created will span the entire CA coast. Use matlab's zoom tool,
% either by typeing 'zoom on' at the command line or selecting the zoom tool from
% the figures toolbar, to select the area of interest.

% Brian Schlining
% 13 Apr 2000

switch nargin
case 0
   fprintf(1, 'No input arguments were specified\n Use as surGrid(sbeFile, <column or dataFile>)\n')
   return
case 1
   fprintf(1, 'Use as surGrid(sbeFile, <column or dataFile>)\n')
   return
end

%*********************
% Load the data files
%*********************
[p f e] = fileparts(sbeFile);
% filename = [p filesep f '.mat'];
filename = [f '.mat'];
if exist(filename) ~= 2
   X   = readPtSurUnderway(sbeFile);              % Load underway data from seabird
   save(filename, 'X');
else
   load(filename);
end

t   = X.time;
lat = X.lat;                     
lon = X.lon;

if isstr(column)                  % If the second input argument is a string it specifes a file of underway data (like pCO2)
   x   = load(column);            % Load the data, (a delimited file saved from an excel sheet)
   tx  = xl2sdn(x(:,1));          % Convert the xl times in the first column to matlab dates
   dat = x(:,2);                  % Second column contains data
   lat = interp1(t, lat, tx);     % Navigate the data base on the seabird data file
   lon = interp1(t, lon, tx);
else
   dat = X.data(:, column);
end


if ~isempty(dat)

   % Marcelino!! Change the collowing line to adust color limits
   CLim = [prctile_(dat, 3) prctile_(dat, 97)];   % Specifies the prctile of data ato use for setting colorLimits

   %***************
   % Grid the data
   %***************
   good = find(~isnan(lat) & ~isnan(lon));
   lat  = lat(good);
   lon  = lon(good);
   dat  = dat(good);
   XI   = [min(lon):(max(lon) - min(lon))./200:max(lon)];
   YI   = [min(lat):(max(lat) - min(lat))./200:max(lat)]';
   ZI   = griddata(lon,lat,dat,XI,YI, 'linear');  % Type options are linear, cubic or nearest

   %***************
   % Create a plot
   %***************
   figure(1); clf
   pcolor(XI,YI,ZI); shading interp; colorbar; hold on;
   load c:\muse\auv\day\ryjo\monitoring_xy; plot(x,y,'k',x,y,'ko')
   % surf(XI,YI,ZI)                % Data Contour
   % set(gca, 'CLim', CLim)
   % shading interp
   % view(2)
   % hold on
   % ZLim = get(gca, 'ZLim');
   % plot3(lon, lat, ones(size(dat))*(max(ZLim)+1), 'k.', 'MarkerSize', .1)  % Plot sample positions
   load CST                      % Load the coastline data
   plot(CST(:,1),CST(:,2),'k');
   plot(lon, lat, 'k.', 'MarkerSize', .1)  % Plot sample positions
   % plot3(CST(:,1), CST(:,2), ones(size(CST(:,1)))*(max(ZLim)+1), 'k') % Plot the coastline
   % hold off
   % colorbar2
   %set(gca,'XLim', [-122.422358651448  -121.766379333516], 'YLim', [36.56187756710288  37.13674116292442])
   axis([-122.5 -121.7 36.5 37]);
   % axis equal
   xlabel('Longitude')
   ylabel('Latitude')
   box on
   
   %[r c] = size(ZI);
   %XI = repmat(XI, r, 1);
   %YI = repmat(YI, 1, c);
   %dat = [XI(:), YI(:) ZI(:)];
   %save temp.txt dat -ascii
end

function y = prctile_(x,p);
%PRCTILE gives the percentiles of the sample in X.
%   Y = PRCTILE(X,P) returns a value that is greater than P percent
%   of the values in X. For example, if P = 50  Y is the median of X.
%
%   P may be either a scalar or a vector. For scalar P, Y is a row
%   vector containing Pth percentile of each column of X. For vector P,
%   the ith row of Y is the P(i) percentile of each column of X.

[prows pcols] = size(p);
if prows ~= 1 & pcols ~= 1
    error('P must be a scalar or a vector.');
end

if any(p > 100) | any(p < 0)
    error('P must take values between 0 and 100');
end

xx = sort(x);
[m,n] = size(x);

if m==1 | n==1
    m = max(m,n);
   if m == 1,
      y = x*ones(length(p),1);
      return;
   end
    n = 1;
    q = 100*(0.5:m - 0.5)./m;
    xx = [min(x); xx(:); max(x)];
else
    q = 100*(0.5:m - 0.5)./m;
    xx = [min(x); xx; max(x)];
end

q = [0 q 100];
y = interp1(q,xx,p);
