cropping netcdf files using geo-coordinates

I want to crop my study area using lat and long mask. the code is below

%Cropping the files
chlFileDirectory = "/Users/gulfcarbon/Desktop/Chlorophyll_MODIS/2021";
chlFilename = dir(fullfile(chlFileDirectory, '*_chl.nc'));
nChlFileDirectory = length(chlFilename);
outputFolder = "/Users/gulfcarbon/Desktop/Chlorophyll_MODIS/Cropped_2021";
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder);
end
% Define the latitude and longitude ranges for cropping
crop_latitude_range = [30, 24];
crop_longitude_range = [-94, -84];
%lat_mask = crop_latitude_range(2);
% ...

for i = 1:nChlFileDirectory
    chl_nc_filename = chlFilename(i).name;
    chl_nc_file_path = fullfile(chlFileDirectory, chl_nc_filename);
    
    % Read chlorophyll data from the current file
    chl_data = ncread(chl_nc_file_path, 'chl_image');
    latitude = ncread(chl_nc_file_path, 'latitude');
    longitude = ncread(chl_nc_file_path, 'longitude');
    
    % Perform logical indexing to crop the data within the specified latitude and longitude ranges
    lat_mask = (latitude >= crop_latitude_range(2)) & (latitude <= crop_latitude_range(1));
    lon_mask = (longitude >= crop_longitude_range(1)) & (longitude <= crop_longitude_range(2));
    
    % Apply the logical masks to obtain the cropped data
    chl_data_cropped = chl_data(lat_mask, lon_mask);
    latitude_cropped = latitude(lat_mask);
    longitude_cropped = longitude(lon_mask);
    
    % Create the output chlorophyll file name with "_cropped" suffix
    [filepath, name, ext] = fileparts(chl_nc_filename);
    chl_nc_filename_cropped = strcat([name, '_cropped', ext]);
    chl_nc_file_path_cropped = fullfile(outputFolder, chl_nc_filename_cropped);
    
    % Write cropped chlorophyll data to a new netCDF file
    nccreate(chl_nc_file_path_cropped, 'chl_image', 'Dimensions', {'row', numel(latitude_cropped), 'col', numel(longitude_cropped)});
    nccreate(chl_nc_file_path_cropped, 'longitude', 'Dimensions', {'row', numel(latitude_cropped), 'col', numel(longitude_cropped)});
    nccreate(chl_nc_file_path_cropped, 'latitude', 'Dimensions', {'row', numel(latitude_cropped), 'col', numel(longitude_cropped)});
    
    ncwrite(chl_nc_file_path_cropped, 'longitude', longitude_cropped);
    ncwrite(chl_nc_file_path_cropped, 'latitude', latitude_cropped);
    ncwrite(chl_nc_file_path_cropped, 'chl_image', chl_data_cropped);
end

THe error is

The logical indices in position 1 contain a true value outside of the array bounds.

I have added the a image of workspace with the question

Answer 1

[Arnab Paul]

MODIS .nc data can be tricky unlike other climate data format

chlFileDirectory = "Your folder";
chlFilename = dir(fullfile(chlFileDirectory, '*_chl.nc'));
nChlFiles = length(chlFilename);  % Rename the variable
outputFolder = "Folder for the output files";
for i = 1:nChlFiles
    chl_nc_filename = chlFilename(i).name;
    chl_nc_file_path = fullfile(chlFileDirectory, chl_nc_filename);
    min_lat = 27;
    max_lat = 31;
    min_lon = -95;
    max_lon = -88;
    % Read chlorophyll data from the current file
    chl_image = ncread(chl_nc_file_path, 'chl_image');
    latitude = ncread(chl_nc_file_path, 'latitude');
    longitude = ncread(chl_nc_file_path, 'longitude');
    
    in_box_idx_lat = find(latitude >= min_lat & latitude <= max_lat) ;
    in_box_idx_lon =find(longitude >= min_lon & longitude <= max_lon);
    % Clip the data to the bounding box
    lat_clip = latitude(in_box_idx_lat);
    lon_clip = longitude(in_box_idx_lon);
    latDiff = abs(latitude(1,1) - latitude(1,2));
    lonDiff = abs(longitude(1,1) - longitude(1,2));
    latLim = max(lat_clip):-latDiff:min(lat_clip); % decsending order
    lonLim = min(lon_clip):lonDiff:max(lon_clip);% ascending order
    
    % Create a meshgrid for desired output points
    [lon_mesh, lat_mesh] = meshgrid(lonLim, latLim);
    chl_image_clipped = griddata(double(latitude), double(longitude), double(chl_image),...
        double(lat_mesh), double(lon_mesh));
    dim_clip = size(chl_image_clipped);
    nrows = dim_clip(1);
    ncols = dim_clip(2);   
    
    % Create the output chlorophyll file name
    [filepath, name, ext] = fileparts(filename(i).name);
    chl_nc_filename = strcat(name, '_clipped_chl.nc');  
    chl_nc_file_path = fullfile(outputFolder, chl_nc_filename);
    nccreate(chl_nc_file_path, 'chl_image', 'Dimensions', {'row', nrows, 'col', ncols});
    nccreate(chl_nc_file_path,  'longitude_clip', 'Dimensions', {'row', nrows, 'col',ncols });
    nccreate(chl_nc_file_path,  'latitude_clip', 'Dimensions', {'row', nrows, 'col', ncols});
    
    ncwrite(chl_nc_file_path, 'longitude_clip', lon_mesh);
    ncwrite(chl_nc_file_path, 'latitude_clip', lat_mesh);
    ncwrite(chl_nc_file_path, 'chl_image', chl_image_clipped);
end

Answer 2

[KSSV]

Why don't you use interp2.

Let X, Y and Z be your original data. Xi, Yi is your coordinates for which you want to crop/ extract the data.

Zi = interp2(X,Y,Z,Xi,Yi) ;

[Arnab Paul]

for i = 1:nChlFileDirectory
    chl_nc_filename = chlFilename(i).name;
    chl_nc_file_path = fullfile(chlFileDirectory, chl_nc_filename);
    
    % Read chlorophyll data from the current file
    chl_data = ncread(chl_nc_file_path, 'chl_image');
    latitude = ncread(chl_nc_file_path, 'latitude');
    longitude = ncread(chl_nc_file_path, 'longitude');
    
    % Define the latitude and longitude ranges for cropping
    crop_latitude_range = [30, 24];
    crop_longitude_range = [-94, -84];
    
    % Create grid for original latitude and longitude data
    [X, Y] = meshgrid(longitude, latitude);
    
    % Create grid for cropped latitude and longitude data
    [Xi, Yi] = meshgrid(crop_longitude_range(1):crop_longitude_range(2), crop_latitude_range(1):-1:crop_latitude_range(2));
    
    % Crop the data using interp2
    chl_data_cropped = interp2(X, Y, chl_data, Xi, Yi, 'linear', NaN);
    

I used this snippet, but the error is

Error using repmat Requested 2748620x2748620 (56288.5GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive.

Error in meshgrid (line 61) xx = repmat(xrow,size(ycol));

[KSSV]

It seems your longitude, latitude are already in meshgrid format. Check their dimensions.

[Arnab Paul]

Both of them are in 1354x2030. That is true for all the images I have in the folder

[KSSV]

That means you need not to use mesh grid.

[Arnab Paul]

crop_latitude_range = [29.75, 27.75];
crop_longitude_range = [-97.30, -87.75];

for i = 1:nChlFileDirectory
    chl_nc_filename = chlFilename(i).name;
    chl_nc_file_path = fullfile(chlFileDirectory, chl_nc_filename);
    
    %% Read chlorophyll data from the current file
    chl_data = ncread(chl_nc_file_path, 'chl_image');
    latitude = ncread(chl_nc_file_path, 'latitude');
    longitude = ncread(chl_nc_file_path, 'longitude');
     % Create grid for cropped latitude and longitude data
    [Xi, Yi] = ndgrid(crop_longitude_range(1):crop_longitude_range(2), crop_latitude_range(1):crop_latitude_range(2));
    % Perform 2D interpolation using griddata
    chl_data_cropped = griddata(longitude, latitude, chl_data, Xi, Yi);
end

Unrecognized function or variable 'nChlFileDirectory'.

I used this but the Xi and Yi gave me a 3x10 matrix with 1 degree gap. not matching with my original chlorophyll image. I have attche the image.