match_roms_mask.m

function output_data = match_roms_mask ( lon, lat, mask, depth, data )
% MATCH ROMS MASK:  make the otps data match the roms land mask.
%
% The data is filled in used the mode of the nearby data rather
% than the mean.  This way the phase doesn't get screwed up as bad.
%
% USAGE:  output_data = match_roms_mask ( lon, lat, mask, data );
%

filled_in_data = data;

%
% Find the ROMS water points (mask == 1) [and water points from grid (depth>21)] which the OTPS grid thought 
% was water (input_constituent == nan).

% SEC change 2/2016 we need to find depth >21 not mask==1 
% nan_inds = find (isnan(data) & mask==1);
%  df = find(lon<-70.75 & lat<43.07); 
 nan_inds = find (isnan(data)==0 & depth>20 & mask==1);

%
% Get the row and column numbers of these points.
[nan_r, nan_c] = indexit_rc ( size(lon), nan_inds );
[num_rows,num_cols] = size(lon);

%
% loop thru each point, trying to fill it in with the nearest non-nan points
num_points = length(nan_inds);

for j = 1:num_points

	%
	% Construct an index mask around this point
	nan_index = nan_inds(j);
	mask_thickness = 1;
	while 1

		nbrows = [nan_r(j)-mask_thickness:nan_r(j)+mask_thickness];
		ind = find((nbrows<1) | (nbrows>num_rows));
		nbrows(ind) = [];

		nbcols = [nan_c(j)-mask_thickness:nan_c(j)+mask_thickness];
		ind = find((nbcols<1) | (nbcols>num_cols));
		nbcols(ind) = [];

		[C,R] = meshgrid(nbcols,nbrows);
		neighbor_inds = (C-1)*num_rows + R;
		%
		% check that none of these neighbor points are land points
% 		valid_ind = find(mask(neighbor_inds) == 1);

% SEC change 2/2016 we need to find depth>21 not mask. We are including land
% in this model 

        % Valid indices are where the depth of the neighbor indices
        % are greater than 21
		valid_ind = find(depth(neighbor_inds)> 21);
        % only keep the data that is greater than 21
		neighbor_inds = neighbor_inds(valid_ind);

		finite_ind = find(isfinite(data(neighbor_inds)) & depth(neighbor_inds)>21);
        
		if any(finite_ind)
			break;
		end

		mask_thickness = mask_thickness + 1;
        
		if ( mask_thickness > 35 )
			error ( 'how can it be this bad?' );
		end
	end

	filled_in_data(nan_index) = my_mode(data(neighbor_inds(finite_ind)));
% 	filled_in_data(nan_index) = mean(data(neighbor_inds(finite_ind)));

end

output_data = filled_in_data;

return

function the_mode = my_mode ( input_data )
mode_data = sort ( input_data );
num_mode_points = length(input_data);
if floor(num_mode_points/2) ~= num_mode_points/2
	mode_index = ceil(num_mode_points/2);
	the_mode = mode_data(mode_index);
else
	mode_index = [num_mode_points/2 num_mode_points/2+1];
	the_mode = mean ( mode_data(mode_index) );
end