Setting section = NaN deletes part of matrix

1 view (last 30 days)
Cecily Tye
Cecily Tye on 20 Jan 2022
Commented: _ on 21 Jan 2022
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;
  5 Comments
Matt J
Matt J on 21 Jan 2022
Perhaps you could show the output of whos() so we can see the sizes of everything in your workspace.

Sign in to comment.

Accepted Answer

_
_ on 21 Jan 2022
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

More Answers (0)

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by