MATLAB Answers

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

4 views (last 30 days)
Karthik M
Karthik M on 14 Sep 2021
Commented: Karthik M on 14 Sep 2021
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
Kindly help, waitng for your valuable inputs
Thanks in advance
regards
  4 Comments
Karthik M
Karthik M on 14 Sep 2021
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

Sign in to comment.

Accepted Answer

KSSV
KSSV on 14 Sep 2021
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])
  5 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!