cropping netcdf files using geo-coordinates

14 次查看(过去 30 天)
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

采纳的回答

Arnab Paul
Arnab Paul 2023-8-30
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

更多回答(1 个)

KSSV
KSSV 2023-8-1
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) ;
  5 个评论
Arnab Paul
Arnab Paul 2023-8-3
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.

请先登录,再进行评论。

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by