Setting section = NaN deletes part of matrix

1 view (last 30 days)
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
Cecily Tye
Cecily Tye on 20 Jan 2022
@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
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

Voss
Voss 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
  2 Comments
Cecily Tye
Cecily Tye on 21 Jan 2022
Omg, knew it would be something dumb like that. Thank you!

Sign in to comment.

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!