Setting section = NaN deletes part of matrix

Hi,

I have a bathymetry dataset and am trying to set a portion of inland lakes that have elevations <0 meters equal to NaN (I'm only interested in the ocean data).

X is a meshgrid matrix of longitude where each column contains the same longitude, and Y is a meshgrid matrix of latitude where each row contains the same latitude. bath_50m is a matrix of the channel with values for depths between 0 and -50 meters and NaNs everywhere else (created with a for loop). X, Y, and bath_50m all have the same size before this error: for some reason, when I run this chunk of code to try and get rid of the lakes, it changes the size of my matrix from 1179x11289 to 7453x11289. I have no idea why it is doing that because I have used this approach many times, and when I try a similar approach for this test matrix there does not seem to be any issues. I double checked that ilon and ilat contain the correct indices I'm looking for. Any idea what the issue is?

% netcdf file of bathymetry 
file = 'santa_barbara_13_mhw_2008.nc';
%ncdisp(file)
lat = ncread(file, 'lat');
lon = ncread(file, 'lon'); 
bath = ncread(file, 'Band1');
lat = reshape(lat,1,[]);

% take a splice of dataset since it is so large 
% let's say approx -120.58116 to -119.46479 for lon 
% and 34.341850 to 34.4846 for lat 
% lon: 
lon1 = -120.58116;
lon2 = -119.46479;
i_lon1 = find(min(abs(lon - lon1)) == abs(lon - lon1));
i_lon2 = find(min(abs(lon - lon2)) == abs(lon - lon2));
ilon = [i_lon1:i_lon2];
% lat: 
lat1 = 34.32;
lat2 = 34.4846;
i_lat1 = find(min(abs(lat - lat1)) == abs(lat - lat1));
i_lat2 = find(min(abs(lat - lat2)) == abs(lat - lat2));
ilat = [i_lat1:i_lat2];

lon_channel = lon(ilon);
lat_channel = lat(ilat);
bath_channel = bath(ilon, ilat);
% make elevations greater than 0 = NaN
bath_channel(bath_channel > 0) = NaN;
bath_channel = bath_channel.'; % transpose 

% visualize data: 
[X,Y] = meshgrid(lon_channel, lat_channel); 

clear lon1 lon2 lat1 lat2 i_lon1 i_lon2 i_lat1 i_lat2 ilon ilat 
%% find places where bathymetry is <= 50 m 
bath_50m = zeros(size(bath_channel));

for i = 1:numel(bath_channel)
    if bath_channel(i) >= -50 
        bath_50m(i) = bath_channel(i);
    else
        bath_50m(i) = NaN; 
    end
end 

lat1 = 34.4169;
lat2 = 34.4435; 
i_lat1 = find(min(abs(Y(:,1) - lat1)) == abs(Y(:,1) - lat1));
i_lat2 = find(min(abs(Y(:,1) - lat2)) == abs(Y(:,1) - lat2));
ilat = [i_lat1:i_lat2];
lon1 = -119.8620;
lon2 = -119.8200; 
i_lon1 = find(min(abs(X(1,:) - lon1)) == abs(X(1,:) - lon1));
i_lon2 = find(min(abs(X(1,:) - lon2)) == abs(X(1,:) - lon2));
ilon = [i_lon1:i_lon2]'; % don't think it's necessary to transpose this, tried it both ways to see if that would help 
bath_50m(ilon, ilat) = NaN; 

% test matrix: 
test_matrix = [1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5];
row = [2:3];
col = [3:4];
test_matrix(row, col) = NaN;

[Matt J]

X, Y, and bath_50m are missing from your post.

[Cecily Tye]

@Matt J sorry, just attached. It is very large

[Matt J]

We still can't run it because you haven't attached santa_barbara_13_mhw_2008 copy.nc

[Cecily Tye]

@Matt J Oh I guess it failed since even the zip is over 5 MB. The data is from this website: https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.ngdc.mgg.dem:603 (netcdf version)

[Matt J]

Perhaps you could show the output of whos() so we can see the sizes of everything in your workspace.

Answer 1

[Voss]

Change this:

bath_50m(ilon, ilat) = NaN;

to this:

bath_50m(ilat, ilon) = NaN;

The reason is because this happened:

 bath_channel = bath_channel.'; % transpose 

[Cecily Tye]

Omg, knew it would be something dumb like that. Thank you!

[Voss]

You are welcome!