How to set level values of a variable by comparing(subtracing) two nc files
1 view (last 30 days)
Show older comments
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
Accepted Answer
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
KSSV
on 14 Sep 2021
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 ;
More Answers (0)
See Also
Categories
Find more on Colormaps in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!