% This is for fixing the caloosahatchee river per Daves request.
% The river is supposed to be shallower near the mouth (6m?) and deeper
% towards the locks. Towards the locks use the bathy that Dave collected.
% The output variable HS is the new bathymetry.
%
%       21 November, 2009

nsmooth = 10;% Number of times to smooth to shallow up the southern portion 
             % of the model grid. The higher the number the shallower it
             % becomes.
crop_d = 10; % Depth to crop the bathymetry at before starting the smoothing
             % to shallow up the bathymetry.
hmin = 3;    % minimum waterdepth
method = 2;  % Set as 1 or 2.
             % Method 1 uses the original bathymetry files to create an
             % entirely new bathymetry. Method 2 uses the smoothed
             % ROMS grid bathymetry made from the original data to generate
             % a new bathymetry. Both methods use the Dave Fugate measured
             % bathy for the northern portion.

% Getting the original grid file name and relavent variables
fname = 'calus_grid1.nc';
rlon = nc_varget(fname,'lon_rho');
rlat = nc_varget(fname,'lat_rho');
mask = nc_varget(fname,'mask_rho');
ho = nc_varget(fname,'h');


% loading in the two sets of bathymetry
B1 = load('swf_bathy_seagrid.mat');   % large bathymetry set
B2 = load('upperCalus.mat');          % Dave measurements (small data set)


if(method == 1);
% remove the data from the large dataset that the small set encompasses
AX = [-81.8621  -81.5827   26.6638   26.7475];   % crop this out of larger dataset
tmp = find((B1.xbathy>AX(1) & B1.xbathy<AX(2)) & ...
    (B1.ybathy>AX(3) & B1.ybathy<AX(4)));
tmp2 = ones(length(B1.xbathy),1);
tmp2(tmp) = 0;
tmp2 = logical(tmp2);

B1.xbathy = B1.xbathy(tmp2);
B1.ybathy = B1.ybathy(tmp2);
B1.zbathy = B1.zbathy(tmp2);
clear tmp tmp2

%Setting a third data set that is  based on the short dataset so the
%interpolation below will fully work.
B3.x = [-81.7471 -81.7617 -81.8085 -81.8091 -81.8259 -81.8394 -81.7545 ...
    -81.7088 -81.7564 -81.8412 -81.6888];
B3.y = [ 26.7185 26.7169 26.6994 26.6995 26.6939 26.6921 26.7226 ...
    26.7271 26.7268 26.6966 26.7318];
B3.z = [4 5.75 5.75 5 4 3 5 6.5 5 3 6.5];

% Combine the three data sets into one.
X = [B1.xbathy;B2.xbathy';B3.x'];
Y = [B1.ybathy;B2.ybathy';B3.y'];
Z = [B1.zbathy;-B2.zbathy';-B3.z'];

% Crop the datasets to speed interpolation. Here, only keeping what is
% needed to interpolate onto the caloosahatchee grid file.
AX = [-82.0482  -81.6716   26.4310   26.7764];
tmp = find((X>AX(1) & X<AX(2)) & (Y>AX(3) & Y<AX(4)));
X = X(tmp);
Y = Y(tmp);
Z = Z(tmp);

% Making sure Z is not smaller than hmin
Z(abs(Z)<hmin) = -hmin;

% Interpolating onto the ROMS Caloosahatchee River grid file.
H = griddata(X,Y,Z,rlon,rlat);
H(~mask) = -hmin;

% Smoothing the grid out and trying to not preserve the deep parts, per
% conversation with Dave. NOTE: not using the normal shapiro2.m filter here
% because I am not trying to conserve volume, but rather shallow up the
% bathymetry. I am trying to preserve the deeper upper portion where the
% actualy point measurements were.
HS = H;
HS(HS>crop_d) = crop_d;       % cropping to help the smoothing shallow up the deeper portions
HS = smooth2d(HS);
HS(~mask) = -hmin;
for i=1:nsmooth
    HS(:,85:end) = H(:,85:end);
    HS = smooth2d(HS,.5,.5);
    HS(~mask) = -hmin;
end
HS(:,86:end) = H(:,86:end);

HS = -HS;  % ROMS wants positive bathymetry numbers.

% End of method 1;
elseif(method == 2)
    
%Setting a third data set that is  based on the short dataset so the
%interpolation below will fully work.
B3.x = [-81.7471 -81.7617 -81.8085 -81.8091 -81.8259 -81.8394 -81.7545 ...
    -81.7088 -81.7564 -81.8412 -81.6888];
B3.y = [ 26.7185 26.7169 26.6994 26.6995 26.6939 26.6921 26.7226 ...
    26.7271 26.7268 26.6966 26.7318];
B3.z = [4 5.75 5.75 5 4 3 5 6.5 5 3 6.5];


% Remove the data from the original ROMS grid that the small set
% encompasses.
AX = [-81.8621  -81.5827   26.6638   26.7475];   % crop this out of larger dataset
tmp = find((rlon>AX(1) & rlon<AX(2)) & (rlat>AX(3) & rlat<AX(4)));
M = min(tmp)-1;

rlon2 = rlon(1:M);
rlat2 = rlat(1:M);
ho2 = ho(1:M);
clear tmp

% Combine the three data sets into one.
X = [rlon2 B2.xbathy B3.x];
Y = [rlat2 B2.ybathy B3.y];
Z = [ho2 B2.zbathy B3.z];

Z(Z<hmin) = hmin;

% Interpolating onto the ROMS Caloosahatchee River grid file.
H = griddata(X,Y,Z,rlon,rlat);
H(~mask) = hmin;

% Smoothing the grid out and trying to not preserve the deep parts, per
% conversation with Dave. NOTE: not using the normal shapiro2.m filter here
% because I am not trying to conserve volume, but rather shallow up the
% bathymetry. I am trying to preserve the deeper upper portion where the
% actualy point measurements were.
HS = H;
HS(HS>crop_d) = crop_d;     % cropping to 15m to help the smoothing shallow up the deeper portions
HS = smooth2d(HS);
HS(~mask) = -hmin;
for i=1:nsmooth
    HS(:,85:end) = H(:,85:end);
    HS = smooth2d(HS,.5,.5);
    HS(~mask) = hmin;
end
HS(:,86:end) = H(:,86:end);

end % if(method)

fname = 'caloos_grid2.nc'
nc_varput(fname,'h',HS);

