How to set level values of a variable by comparing(subtracing) two nc files

Hi Freinds

I need to spot seasonal trends of chlor_a nc data for last 10 years for which I considered 2 year(2015 and 2016) nc files (Feb)month at the moment

I find error while subtracting both files to compare and identify max persistence (chlor_a conc.) hotspot of particular month during successive year . After executing need to save both files in nc format, [+ve] value means max chlor_a conc. and [-ve] value means less conc.

I have extracted the chlor_a data and have plot with respective x and y coordinates,

I used below code for both nc files

ncfile_1 = 'A20150322015059.L3m_MO_CHL_chlor_a_4km.nc';
lat_1 = ncread(ncfile_1,'lat') ;
lon_1 = ncread(ncfile_1,'lon') ;
chlor_a_1 = ncread(ncfile_1,'chlor_a');
[X_1,Y_1] = meshgrid(lon_1,lat_1) ;
nn = 1000;
xi_1 = linspace(30,100,nn) ;
yi_1 = linspace(0,30,nn) ;
[Xi_1,Yi_1] = meshgrid(xi_1,yi_1);
iwant_1 = interp2(X_1,Y_1,chlor_a_1',Xi_1,Yi_1)' ;
ind_1 = find(iwant_1>10 | iwant_1<0);
%ind = find(iwant_2 > 0 & iwant_1<10);
iwant_1(ind_1) = NaN;
%[min_val_1,min_ind_1] = min(iwant_1',[],'all','linear');
[min_val_1,min_ind_1] = min(reshape(iwant_1',[],1));
[r1_1,c1_1]=ind2sub(nn,min_ind_1);
x_min_1 = xi_1(c1_1);
y_min_1 = yi_1(r1_1);
%[max_val_1,max_ind_1] = max(iwant_1',[],'all','linear');
[max_val_1,max_ind_1] = max(reshape(iwant_1',[],1));
[r2_1,c2_1]=ind2sub(nn,max_ind_1);
x_max_1 = xi_1(c2_1);
y_max_1 = yi_1(r2_1);
figure(1)
pcolor(xi_1,yi_1,iwant_1'); shading interp;
c = colorbar;
cmap = jet(255);
cmap(:,2) = 0;
colormap(cmap)
caxis([0, 5])
hold on
plot(x_min_1,y_min_1,'k +', 'MarkerSize', 10);
plot(x_max_1,y_max_1,'m +', 'MarkerSize', 10);
hold off
legend('','min','max');

I need to get the plot as such after subtraction, please find the attached image

similary need to work on to compare weekly data

data link

https://drive.google.com/file/d/1hAeVz54bKLA1dDyk5uxMyjJIppG-XS77/view?usp=sharing

https://drive.google.com/file/d/1XIcbaEbqLQOPgi0Ooq72J2elvZwfQnoO/view?usp=sharing

Kindly help, waitng for your valuable inputs

Thanks in advance

regards

[KSSV]

I find error while subtracting both files to compare

What error you found? Show us the full error. Which line the error was?

[Karthik M]

Dear @KSSV

Thank you in advance

After indexed(set level value of a variable) I used the above code but couldn't able to subtract as mat or nc file. Actually I need get the final ouput as that sample image attached after subtracting both nc files (idexing)

I have used this below code for subtraction the 2 nc files

ncfile_1 = 'A20150322015059.L3m_MO_CHL_chlor_a_4km.nc';
lat_1 = ncread(ncfile_1,'lat') ;
lon_1 = ncread(ncfile_1,'lon') ;
chlor_a_1 = ncread(ncfile_1,'chlor_a');
[X_1,Y_1] = meshgrid(lon_1,lat_1) ;
xi_1 = linspace(30,100,1000) ;
yi_1 = linspace(0,30,1000) ;
[Xi_1,Yi_1] = meshgrid(xi_1,yi_1);
iwant_1 = interp2(X_1,Y_1,chlor_a_1',Xi_1,Yi_1)' ;

figure(1)
pcolor(xi_1,yi_1,iwant_1'); shading interp;
c = colorbar;
cmap = jet(255);
cmap(:,2) = 0;
colormap(cmap)
caxis([0, 5])

ncfile_2 = 'A20160322016060.L3m_MO_CHL_chlor_a_4km.nc';
lat_2 = ncread(ncfile_2,'lat') ;
lon_2 = ncread(ncfile_2,'lon') ;
chlor_a_2 = ncread(ncfile_2,'chlor_a');
[X_2,Y_2] = meshgrid(lon_2,lat_2) ;
xi_2 = linspace(30,100,1000) ;
yi_2 = linspace(0,30,1000) ;
[Xi_2,Yi_2] = meshgrid(xi_2,yi_2);
iwant_2 = interp2(X_2,Y_2,chlor_a_2',Xi_2,Yi_2)' ;

figure(2)
pcolor(xi_2,yi_2,iwant_2'); shading interp;
c = colorbar;
cmap = jet(255);
cmap(:,2) = 0;
colormap(cmap)
caxis([0, 5])

figure(3)
diff = iwant_2-iwant_1;
%M = max(diff, [], 'all');
%M
pcolor(diff'); shading interp;
grid on;
imshow(diff);
imgdiff( 4, 3, : ) = [0 0 0];
mask = sum(imgdiff, 3 ) > 0;
 h = imagesc(imgdiff);
h.AlphaData = mask;
%axis([0 30,30 100]);
c = colorbar;
cmap = jet(255);
cmap(:,2) = 0;
colormap(cmap)
caxis([0, 5])

I want to subtract as mat file not as an image, found error as (background image) further need to mask to remove them..

kindly help

regards

[KSSV]

You have not specified what error you are getting.

[Karthik M]

Dear @KSSV

without setting the index found the error image as attached. I want to set level value of variable with >10 and <0 = to ignore the unwanted values with the line in my code [ind = find(iwant_2>10 | iwant_2<0)]; as shown

While subtracting as 2 mat nc file_1 and ncfile_2 (with index limit). I am trying to obtain the sample image attached from the below code

ncfile_1 = 'A20150322015059.L3m_MO_CHL_chlor_a_4km.nc';

lat_1 = ncread(ncfile_1,'lat') ;

lon_1 = ncread(ncfile_1,'lon') ;

chlor_a_1 = ncread(ncfile_1,'chlor_a');

[X_1,Y_1] = meshgrid(lon_1,lat_1) ;

nn = 1000;

xi_1 = linspace(30,100,nn) ;

yi_1 = linspace(0,30,nn) ;

[Xi_1,Yi_1] = meshgrid(xi_1,yi_1);

iwant_1 = interp2(X_1,Y_1,chlor_a_1',Xi_1,Yi_1)' ;

ind_1 = find(iwant_1>10 | iwant_1<0);

%ind = find(iwant_2 > 0 & iwant_1<10);

iwant_1(ind_1) = NaN;

%[min_val_1,min_ind_1] = min(iwant_1',[],'all','linear');

[min_val_1,min_ind_1] = min(reshape(iwant_1',[],1));

[r1_1,c1_1]=ind2sub(nn,min_ind_1);

x_min_1 = xi_1(c1_1);

y_min_1 = yi_1(r1_1);

%[max_val_1,max_ind_1] = max(iwant_1',[],'all','linear');

[max_val_1,max_ind_1] = max(reshape(iwant_1',[],1));

[r2_1,c2_1]=ind2sub(nn,max_ind_1);

x_max_1 = xi_1(c2_1);

y_max_1 = yi_1(r2_1);

figure(1)

pcolor(xi_1,yi_1,iwant_1'); shading interp;

c = colorbar;

cmap = jet(255);

cmap(:,2) = 0;

colormap(cmap)

caxis([0, 5])

hold on

plot(x_min_1,y_min_1,'k +', 'MarkerSize', 10);

plot(x_max_1,y_max_1,'m +', 'MarkerSize', 10);

hold off

legend('','min','max');

Image = getframe(gcf);

imwrite(Image.cdata, 'feb2015.tif');

--

ncfile_2 = 'A20160322016060.L3m_MO_CHL_chlor_a_4km.nc';

lat_2 = ncread(ncfile_2,'lat') ;

lon_2 = ncread(ncfile_2,'lon') ;

chlor_a_2 = ncread(ncfile_2,'chlor_a');

[X_2,Y_2] = meshgrid(lon_2,lat_2) ;

nn = 1000;

xi_2 = linspace(30,100,nn) ;

yi_2 = linspace(0,30,nn) ;

[Xi_2,Yi_2] = meshgrid(xi_2,yi_2);

iwant_2 = interp2(X_2,Y_2,chlor_a_2',Xi_2,Yi_2)' ;

ind = find(iwant_2>10 | iwant_2<0);

%ind = find(iwant_2 > 0 & iwant_2<10);

iwant_2(ind) = NaN;

%[min_val,min_ind] = min(iwant_2',[],'all','linear');

[min_val,min_ind] = min(reshape(iwant_2',[],1));

[r1,c1]=ind2sub(nn,min_ind);

x_min = xi_2(c1);

y_min = yi_2(r1);

%[max_val,max_ind] = max(iwant_2',[],'all','linear');

[max_val,max_ind] = max(reshape(iwant_2',[],1));

[r2,c2]=ind2sub(nn,max_ind);

x_max = xi_2(c2);

y_max = yi_2(r2);

figure(2)

pcolor(xi_2,yi_2,iwant_2'); shading interp;

c = colorbar;

cmap = jet(255);

cmap(:,2) = 0;

colormap(cmap)

caxis([0, 5])

hold on

plot(x_min,y_min,'k +', 'MarkerSize', 10);

plot(x_max,y_max,'m +', 'MarkerSize', 10);

hold off

legend('','min','max');

I have tried my level best for subtracting the 2 files after indexing. I don't know how to execute

kindly help, need to execute file as sample image

thank you

regards

Answer 1

[KSSV]

Your interpolation grid coordinates are completely different comapred to original grid. Read about interp2.

ncfile_1 = 'A20150322015059.L3m_MO_CHL_chlor_a_4km.nc';
lat_1 = ncread(ncfile_1,'lat') ;
lon_1 = ncread(ncfile_1,'lon') ;
chlor_a_1 = ncread(ncfile_1,'chlor_a');
[X_1,Y_1] = meshgrid(lon_1,lat_1) ;

xi_1 = linspace(min(lon_1),max(lon_1),1000) ;   % Changed here 
yi_1 = linspace(min(lat_1),max(lat_1),1000) ;   % changed here 

[Xi_1,Yi_1] = meshgrid(xi_1,yi_1);
iwant_1 = interp2(X_1,Y_1,chlor_a_1',Xi_1,Yi_1)' ;

figure(1)
pcolor(xi_1,yi_1,iwant_1'); shading interp;
c = colorbar;
cmap = jet(255);
cmap(:,2) = 0;
colormap(cmap)
caxis([0, 5])

ncfile_2 = 'A20160322016060.L3m_MO_CHL_chlor_a_4km.nc';
lat_2 = ncread(ncfile_2,'lat') ;
lon_2 = ncread(ncfile_2,'lon') ;
chlor_a_2 = ncread(ncfile_2,'chlor_a');
[X_2,Y_2] = meshgrid(lon_2,lat_2) ;

xi_2 = linspace(min(lon_2),max(lon_2),1000) ;   % Changed here 
yi_2 = linspace(min(lat_2),max(lat_2),1000) ;   % Changed here 

[Xi_2,Yi_2] = meshgrid(xi_2,yi_2);
iwant_2 = interp2(X_2,Y_2,chlor_a_2',Xi_2,Yi_2)' ;

figure(2)
pcolor(xi_2,yi_2,iwant_2'); shading interp;
c = colorbar;
cmap = jet(255);
cmap(:,2) = 0;
colormap(cmap)
caxis([0, 5])

figure(3)
diff = iwant_2-iwant_1;
%M = max(diff, [], 'all');
%M
pcolor(diff'); shading interp;
grid on;
imshow(diff);
imgdiff( 4, 3, : ) = [0 0 0];
mask = sum(imgdiff, 3 ) > 0;
 h = imagesc(imgdiff);
h.AlphaData = mask;
%axis([0 30,30 100]);
c = colorbar;
cmap = jet(255);
cmap(:,2) = 0;
colormap(cmap)
caxis([0, 5])

[Karthik M]

Dear @KSSV

without setting the index found the error image as attached

I want to set level value of variable with >10 and <0 = to ignore the values

with the line in my code [ind = find(iwant_2>10 | iwant_2<0)]; as shown

While subtracting as 2 mat nc file_1 and ncfile_2 (with index limit)

Need to obtain the sample image attached from the code below

ncfile_1 = 'A20150322015059.L3m_MO_CHL_chlor_a_4km.nc';

lat_1 = ncread(ncfile_1,'lat') ;

lon_1 = ncread(ncfile_1,'lon') ;

chlor_a_1 = ncread(ncfile_1,'chlor_a');

[X_1,Y_1] = meshgrid(lon_1,lat_1) ;

nn = 1000;

xi_1 = linspace(30,100,nn) ;

yi_1 = linspace(0,30,nn) ;

[Xi_1,Yi_1] = meshgrid(xi_1,yi_1);

iwant_1 = interp2(X_1,Y_1,chlor_a_1',Xi_1,Yi_1)' ;

ind_1 = find(iwant_1>10 | iwant_1<0);

%ind = find(iwant_2 > 0 & iwant_1<10);

iwant_1(ind_1) = NaN;

%[min_val_1,min_ind_1] = min(iwant_1',[],'all','linear');

[min_val_1,min_ind_1] = min(reshape(iwant_1',[],1));

[r1_1,c1_1]=ind2sub(nn,min_ind_1);

x_min_1 = xi_1(c1_1);

y_min_1 = yi_1(r1_1);

%[max_val_1,max_ind_1] = max(iwant_1',[],'all','linear');

[max_val_1,max_ind_1] = max(reshape(iwant_1',[],1));

[r2_1,c2_1]=ind2sub(nn,max_ind_1);

x_max_1 = xi_1(c2_1);

y_max_1 = yi_1(r2_1);

figure(1)

pcolor(xi_1,yi_1,iwant_1'); shading interp;

c = colorbar;

cmap = jet(255);

cmap(:,2) = 0;

colormap(cmap)

caxis([0, 5])

hold on

plot(x_min_1,y_min_1,'k +', 'MarkerSize', 10);

plot(x_max_1,y_max_1,'m +', 'MarkerSize', 10);

hold off

legend('','min','max');

Image = getframe(gcf);

imwrite(Image.cdata, 'feb2015.tif');

--

ncfile_2 = 'A20160322016060.L3m_MO_CHL_chlor_a_4km.nc';

lat_2 = ncread(ncfile_2,'lat') ;

lon_2 = ncread(ncfile_2,'lon') ;

chlor_a_2 = ncread(ncfile_2,'chlor_a');

[X_2,Y_2] = meshgrid(lon_2,lat_2) ;

nn = 1000;

xi_2 = linspace(30,100,nn) ;

yi_2 = linspace(0,30,nn) ;

[Xi_2,Yi_2] = meshgrid(xi_2,yi_2);

iwant_2 = interp2(X_2,Y_2,chlor_a_2',Xi_2,Yi_2)' ;

ind = find(iwant_2>10 | iwant_2<0);

%ind = find(iwant_2 > 0 & iwant_2<10);

iwant_2(ind) = NaN;

%[min_val,min_ind] = min(iwant_2',[],'all','linear');

[min_val,min_ind] = min(reshape(iwant_2',[],1));

[r1,c1]=ind2sub(nn,min_ind);

x_min = xi_2(c1);

y_min = yi_2(r1);

%[max_val,max_ind] = max(iwant_2',[],'all','linear');

[max_val,max_ind] = max(reshape(iwant_2',[],1));

[r2,c2]=ind2sub(nn,max_ind);

x_max = xi_2(c2);

y_max = yi_2(r2);

figure(2)

pcolor(xi_2,yi_2,iwant_2'); shading interp;

c = colorbar;

cmap = jet(255);

cmap(:,2) = 0;

colormap(cmap)

caxis([0, 5])

hold on

plot(x_min,y_min,'k +', 'MarkerSize', 10);

plot(x_max,y_max,'m +', 'MarkerSize', 10);

hold off

legend('','min','max');

kindly help, need to execute file as sample image

thank you

regards

[KSSV]

Replace unwanted values with NaN's by indexing.

[Karthik M]

Dear @KSSV

I have tried my level best for subtracting the 2 files after indexing. I don't know how to execute, kindly help

regards

[KSSV]

If A is your data matrix; to ingnore values which are <0 and >10; replace them with NaN's using:

idx = A < 0 | A > 10 ;   % this will give logical indices 
A(idx) = NaN ; 

[Karthik M]

Thank you @KSSV for your valuable time