Cannot write netcdf file

filename = "/Users/gulfcarbon2/Desktop/AD-ATCOR/output/AQUA_MODIS.20190816T194000.L2._Rrs.nc";
info = ncinfo(filename);
dim = cell2mat({info.Variables.Dimensions});
ncol = dim(1).Length;
nrow = dim(2).Length;
% Read the variables from the NetCDF file
Rrs412 = ncread(filename, 'Rrs412');
Rrs443 = ncread(filename, 'Rrs443');
Rrs469 = ncread(filename, 'Rrs469');
Rrs488 = ncread(filename, 'Rrs488');
Rrs531 = ncread(filename, 'Rrs531');
Rrs547 = ncread(filename, 'Rrs547');
Rrs555 = ncread(filename, 'Rrs555');
Rrs645 = ncread(filename, 'Rrs645');
Rrs667 = ncread(filename, 'Rrs667');
Rrs678 = ncread(filename, 'Rrs678');
Rrs748 = ncread(filename, 'Rrs748');
Rrs859 = ncread(filename, 'Rrs859');
Rrs869 = ncread(filename, 'Rrs869');
latitude = ncread(filename, 'latitude');
longitude = ncread(filename, 'longitude');

%Threshold 
rho = log(Rrs555 ./ Rrs667);
% Calculate X as log[(max(Rrs443, Rrs489) / Rrs555)] for OC3M algorithm
X = log(max(Rrs443, Rrs488) ./ Rrs555);

% Initialize an array for the resulting pixels
result_image = zeros(size(longitude));

% Apply the condition and equations
a0 = 0.283;
a1 = -2.753;
a2 = 1.457;
a3 = 0.659;
a4 = -1.403;
% Apply the condition and equations
threshold = 0.75;

% Iterate over each pixel
for i = 1:numel(rho)
    if rho(i) < threshold
        % Apply the equation pixel = 0.3894 * rho + 2.9642 for rho < threshold
        result_image(i) = 0.3894 * rho(i) + 2.9642;
    else
        % Apply the equation y = exp(a0 + a1 * X + a2 * X^2 + a3 * X^3 + a4 * X^4) for rho >= threshold
        result_image(i) = exp(a0 + a1 * X(i) + a2 * X(i)^2 + a3 * X(i)^3 + a4 * X(i)^4);
    end
end

% % Plot the result on a map
% Ad_fig = figure;
% m_proj('Mercator','longitude',([round(min(longitude(:,1)))-0.5 round(max(longitude(:,1)))+0.5]),'latitude',([round(min(latitude(1,:)))-0.5 round(max(latitude(1,:)))+0.5]));
% m_pcolor(longitude, latitude, result_image);
% shading interp;
% m_gshhs_f('patch',[.5 .5 .5],'edgecolor','none');
% m_grid('linewi',2,'tickdir','out','box','fancy','fontsize',16);
% m_ruler([.7 1],.9,'tickdir','out','ticklen',[.007 .007],'fontsize',12);
% h=colorbar;
% % set(get(h,'ylabel'),'String',sprintf(result_image),'fontsize',16);
% colormap(m_colmap('jet'));
% caxis([0 0.03]);
% brighten(.2)
% m_coast('color', 'k');
% title('Result Image');

% Create a new NetCDF file for the result
nc_filename = [filename(1:end-6), 'AD_Rrs.nc'];

% Write the result data to the NetCDF file
nccreate(nc_filename, 'result_image', 'Dimensions', {'r', nrow, 'c', ncol}, 'Format', 'classic');
ncwrite(nc_filename, 'result_image', result_image);
nccreate(nc_filename, 'latitude', 'Dimensions', {'r', nrow, 'c', ncol}, 'Format', 'classic');
ncwrite(nc_filename, 'latitude', latitude);
nccreate(nc_filename, 'longitude', 'Dimensions', {'r', nrow, 'c', ncol}, 'Format', 'classic');
ncwrite(nc_filename, 'longitude', longitude);

So far I have done to get a plot. Now I want to save it as a netcdf file. I can also save the the .nc file but while opening it in a software such as seadas the longitude and latitude data is missing. error coming such as

Error using nccreate
'result_image' exists in file. Can not modify existing variables.

Error in AD_Arnab (line 75)
nccreate(nc_filename, 'result_image', 'Dimensions', {'r', nrow, 'c', ncol}, 'Format', 'classic');

any help will be appriciated.

Answer 1

[dpb]

filename = "/Users/gulfcarbon2/Desktop/AD-ATCOR/output/AQUA_MODIS.20190816T194000.L2._Rrs.nc"

filename = "/Users/gulfcarbon2/Desktop/AD-ATCOR/output/AQUA_MODIS.20190816T194000.L2._Rrs.nc"

nc_filename = [filename(1:end-6), 'AD_Rrs.nc']

nc_filename = "AD_Rrs.nc"

nrow=10; ncol=10;
result_image=rand(nrow,ncol);
nccreate(nc_filename, 'result_image', 'Dimensions', {'r', nrow, 'c', ncol}, 'Format', 'classic');
ncwrite(nc_filename, 'result_image', result_image);
ncdisp(nc_filename)

Source:
           /users/mss.system.BSKLtc/AD_Rrs.nc
Format:
           classic
Dimensions:
           r = 10
           c = 10
Variables:
    result_image
           Size:       10x10
           Dimensions: r,c
           Datatype:   double

Basic code/function seems to work here; we can't duplicate your identical case without the file, but that doesn't seem to be the issue.

NOTA BENE however, that

nc_filename = [filename(1:end-6), 'AD_Rrs.nc']

nc_filename = "AD_Rrs.nc"

isn't doing what you probably think it is -- referencing the substring out of a string variable requires using the curlies {} just as with a cell to convert the string back to char and subscript within it -- parenetheses refer to elements of an array and the dimension of filename is an array of 1, so (1:end-6) --> (1:1-6) --> (1:-5) which is a null operation.

nc_filename = [filename{:}(1:end-6), 'AD_Rrs.nc']

nc_filename = '/Users/gulfcarbon2/Desktop/AD-ATCOR/output/AQUA_MODIS.20190816T194000.L2._AD_Rrs.nc'

Better yet than counting and coding explicit numbers into a string expression, use

[p,n]=fileparts(filename{:})

p = '/Users/gulfcarbon2/Desktop/AD-ATCOR/output'

n = 'AQUA_MODIS.20190816T194000.L2._Rrs'

nc_filename=fullfile(p,[n 'AD_RRs.nc'])

nc_filename = '/Users/gulfcarbon2/Desktop/AD-ATCOR/output/AQUA_MODIS.20190816T194000.L2._RrsAD_RRs.nc'

As another side note regarding using MATLAB, you don't need to write loops for such operations as done above.

threshold = 0.75;
% Iterate over each pixel

with vector operations in MATLAB is

ixLo=(rho<threshold)        % logical condition addressing array
% Apply the equation pixel = 0.3894 * rho + 2.9642 for rho < threshold
result_image=zeros(size(rho));      % preallocate output
result_image(ixLo) = 0.3894*rho(ixLo) + 2.9642;
% Apply the equation y = exp(a0 + a1 * X + a2 * X^2 + a3 * X^3 + a4 * X^4) for rho >= threshold
result_image(~ixLo) = exp(a0 + a1*X(~ixLO) + a2*X(i~ixLO)^2 + a3*X(~ixLO)^3 + a4*X(~ixLO)^4);

Alternatively, since all elements will be written, you can save one step in the prealocation by simply assigning all to the first condition, then overwriting the alternate...

ixLo=(rho<threshold)        % logical condition addressing array
% Apply the equation pixel = 0.3894 * rho + 2.9642 for rho < threshold
result_image= 0.3894*rho+2.9642;        % assign all first
% Apply the equation y = exp(a0 + a1 * X + a2 * X^2 + a3 * X^3 + a4 * X^4) for rho >= threshold
result_image(~ixLo) = exp(a0 + a1*X(~ixLO) + a2*X(i~ixLO)^2 + a3*X(~ixLO)^3 + a4*X(~ixLO)^4);

will produce the same end result.

One additional MATLAB idiom would be to not write polynomial coefficients as separate variables sequentially numbered but use an array...

a=polyfit(X,log(Y),4);                          % compute the coefficients, a
result_image(~ixLo)=exp(polyval(a,X(~ixLO));    % evaluate it where needed...

NOTA BENE SECUNDUS: The order of coefficients returned by polyfit and expected by polyval is highest order first; ergo if the values are coming from somewhere else as the separate ai, then

a=[a4 a3 a2 a1 a0];

[Arnab Paul]

This looks okay and it worked. What I am looking for is the Ad image to come with latitude and longitude. Which is not working in my code. Do you need specific information of the code and image? So it will become easier for you to understand?

[dpb]

Error using nccreate
'result_image' exists in file. Can not modify existing variables.

Error in AD_Arnab (line 75)
nccreate(nc_filename, 'result_image', 'Dimensions', {'r', nrow, 'c', ncol}, 'Format', 'classic');

says the variable 'result_image' was already in the file; the conclusion must be that the file your code created already existed and had previously been written; quite possible to have occurred during code development/debugging.

My guess is the issue with the file name generation above means the file you think you're creating isn't the actual file and that THAT file does already have content in it.

You'll need to search for all the .nc files that are in the working path and make sure there isn't a culprit already there. If the intent is to create a new file; then you can add code to test if the new filename created already exists as a file before continuting to try to nccreate() into it.

The code above IN ISOLATION will not cause the problem...