I have a problem with my script. I would like to extract pH data from the coast of Cameroon or the Gulf of Guinea. Can I have some ideas on what is bugging my script and a cor
    4 Ansichten (letzte 30 Tage)
  
       Ältere Kommentare anzeigen
    
    Victor Ebolo Nkongo
 am 26 Nov. 2023
  
    
    
    
    
    Kommentiert: Mathieu NOE
      
 am 27 Nov. 2023
            clear all close all clc
disp('------------------------------------------------------') disp('Step 1 : Choix du fichier à traiter') [FileName,PathName,~] = uigetfile; disp('...... OK') disp('------------------------------------------------------')
disp('------------------------------------------------------') disp('Step 2 : Affichage des métadonnées netcdf') ncdisp([PathName,FileName]); disp('...... OK') disp('------------------------------------------------------')
disp('------------------------------------------------------') disp('Step 3 : Extraction des variables à Kribi (3°N; 9.5°E)') lon=ncread([PathName,FileName],'LONGITUDE'); % longitude lat=ncread([PathName,FileName],'LATITUDE'); % Latitude
[~,ilatK]=min(abs(lat-3.4)); % Indice latitude correspondant à Kribi [~,ilonK]=min(abs(lon-6.33)); % Indice longitude correspondant à Kribi
if (~isempty(ilatK)&&~isempty(ilonK)) pH=double(ncread([PathName,FileName],'PHPH',[ilonK,ilatK,1],[1 1 Inf],[1 1 1])); hsK=NaN(1,size(pH,2)); % Initialisation du vecteur de sorties des pH for i=1:size(pH,2) pHK(i)=pH(:,:,i); % Compilation de la série temporelle end
    PH=ncread([PathName,FileName],'PHPH_QC',[ilonK,ilatK,1],[1 1 Inf],[1 1 1]);
    PHK=NaN(1,size(PH,3)); % Initialisation du vecteur de sorties chl
    for i=1:size(PH,3)
        PHK(i)=chl(:,:,i); % Compilation de la série temporelle
    end
end
disp('Step 4 : Time vector') t1=double(ncread([PathName,FileName],'TIME')); % Extraction du vecteur temps t=NaN(1,length(t1)); % Initialisation du vecteur de sorties du temps for i=1:length(t1) t(i)=datenum(1900,1,1,0,0,0)+(t1(i)/24); % Conversion en jours juliens end disp('...... OK')
hold on plot(t,pHK,'.-k','Linewidth',1,'Markersize',10) % Hauteur au large plot(t,PHK,'.-r','Linewidth',1,'Markersize',10) % Hauteur au déferlement hold off config_xdate(t,1,5,'pH',0) ylabel('pH (m)','fontweight','bold','fontsize',15) legend('SEA WATER','quality flag')
0 Kommentare
Akzeptierte Antwort
  Mathieu NOE
      
 am 27 Nov. 2023
        hello 
the provided nc file contains a lot of NaN data so don't be surprised here to get most of the time NaN as a  result to your ncread calls 
otherwise I correected a few bugs (original code is commented and new code is the line right after
clear all 
close all
clc
disp('------------------------------------------------------')
disp('Step 1 : Choix du fichier à traiter')
[FileName,PathName,~] = uigetfile('*.nc');
disp('...... OK')
disp('------------------------------------------------------')
disp('------------------------------------------------------')
disp('Step 2 : Affichage des métadonnées netcdf') 
ncdisp([PathName,FileName]);
disp('...... OK')
disp('------------------------------------------------------')
disp('------------------------------------------------------')
disp('Step 3 : Extraction des variables à Kribi (3°N; 9.5°E)')
lon=ncread([PathName,FileName],'LONGITUDE');% longitude 
lat=ncread([PathName,FileName],'LATITUDE');% Latitude
[~,ilatK]=min(abs(lat-3.4));% Indice latitude correspondant à Kribi 
[~,ilonK]=min(abs(lon-6.33));% Indice longitude correspondant à Kribi
if (~isempty(ilatK)&&~isempty(ilonK))
%     pH=double(ncread([PathName,FileName],'PHPH',[ilonK,ilatK,1],[1 1 Inf],[1 1 1]));
    pH=ncread([PathName,FileName],'PHPH',[ilonK ilatK],[1 1]);
    hsK=NaN(1,size(pH,2)); % Initialisation du vecteur de sorties des pH 
    for i=1:size(pH,2)
        pHK(i)=pH(:,:,i); % Compilation de la série temporelle
    end
%     PH=ncread([PathName,FileName],'PHPH_QC',[ilonK,ilatK,1],[1 1 Inf],[1 1 1]);
    PH=ncread([PathName,FileName],'PHPH_QC',[ilonK ilatK],[1 1]);
    PHK=NaN(1,size(PH,3)); % Initialisation du vecteur de sorties chl
    for i=1:size(PH,3)
        PHK(i)=chl(:,:,i); % Compilation de la série temporelle
    end
end
disp('Step 4 : Time vector')
t1=double(ncread([PathName,FileName],'TIME')); 
% Extraction du vecteur temps 
t=NaN(1,length(t1)); % Initialisation du vecteur de sorties du temps 
for i=1:length(t1)
    t(i)=datenum(1900,1,1,0,0,0)+(t1(i)/24); % Conversion en jours juliens
end
disp('...... OK')
hold on 
plot(t,pHK,'.-k','Linewidth',1,'Markersize',10) % Hauteur au large
plot(t,PHK,'.-r','Linewidth',1,'Markersize',10) % Hauteur au déferlement 
hold off 
config_xdate(t,1,5,'pH',0)
ylabel('pH (m)','fontweight','bold','fontsize',15) 
legend('SEA WATER','quality flag')
3 Kommentare
  Mathieu NOE
      
 am 27 Nov. 2023
				ok , I fixed minors bugs , like 
these 2 lines where missing : 
[~,ilatK]=min(abs(lat-3.4));% Indice latitude correspondant à Kribi 
[~,ilonK]=min(abs(lon-6.33));% Indice longitude correspondant à Kribi
also you can do this
% Extraction du vecteur temps 
t=NaN(1,length(t1)); % Initialisation du vecteur de sorties du temps 
for i=1:length(t1)
    t(i)=datenum(1950,1,1,0,0,0)+(t1(i)/24); % Conversion en jours juliens
end
 in one line without a for loop (dtanum operates on vectors) 
t=datenum(1950,1,1,0,0,0)+(t1/24); % Conversion en jours juliens
full code - but same remark as above, with the provided nc fill all data at the requested lon / lat coordinates are NaNs
the last portion of code I cannot really test it as sss is NaN , so this line of code will fail 
surf(lo,la,sss(:,:,124))
what is the intention as sss is supposed to be a single value ? 
sss=ncread([PathName,FileName],'PSAL',[ilonK ilatK],[1 1]);   will give value PSAL at ilonK ilatK coordinates = 1 value , not the full map
clear all 
close all
clc
disp('------------------------------------------------------')
disp('Step 1 : Choix du fichier à traiter')
[FileName,PathName,~] = uigetfile('*.nc');
disp('...... OK')
disp('------------------------------------------------------')
disp('------------------------------------------------------')
disp('Step 2 : Affichage des métadonnées netcdf') 
ncdisp([PathName,FileName]);
disp('...... OK')
disp('------------------------------------------------------')
disp('------------------------------------------------------')
disp('Step 3 : Extraction des variables à Kribi (3°N; 9.5°E)')
lon=ncread([PathName,FileName],'LONGITUDE');% longitude 
lat=ncread([PathName,FileName],'LATITUDE');% Latitude
[~,ilatK]=min(abs(lat-3.4));% Indice latitude correspondant à Kribi 
[~,ilonK]=min(abs(lon-6.33));% Indice longitude correspondant à Kribi
if (~isempty(ilatK)&&~isempty(ilonK))
    sss=ncread([PathName,FileName],'PSAL',[ilonK ilatK],[1 1]);
    SSSK=NaN(1,size(sss,2)); % Initialisation du vecteur de sorties des pH 
    for i=1:size(sss,2)
        SSSK(i)=sss(:,:,i); % Compilation de la série temporelle
    end
    sal=ncread([PathName,FileName],'PSAL_QC',[ilonK ilatK],[1 1]);
    salK=NaN(1,size(sal,3)); % Initialisation du vecteur de sorties chl
    for i=1:size(sal,3)
        salK(i)=sal(:,:,i); % Compilation de la série temporelle
    end
end
disp('Step 4 : Time vector')
t1=ncread([PathName,FileName],'TIME'); 
% Extraction du vecteur temps 
t=datenum(1950,1,1,0,0,0)+(t1/24); % Conversion en jours juliens
% the beginning of the sequel
[la,lo]=meshgrid(lat,lon);
figure;
title('Temperature')
hold on
surf(lo,la,sss(:,:,124))
view(2)
colorbar
shading interp
plot(9.87,2.95,'.k','MarkerSize',50) % Kribi
plot3(9.25,3.25,1000,'.k','MarkerSize',50) % Point extraction
hold off
  Mathieu NOE
      
 am 27 Nov. 2023
				looking at data information provided by this code
ncdisp([PathName,FileName]);
one extract , for example : 
    PSAL                    
           Size:       2174x903
           Dimensions: DEPTH,TIME
           Datatype:   int32
           Attributes:
                       _FillValue          = -2147483647
                       standard_name       = 'sea_water_practical_salinity'
                       long_name           = 'Practical salinity'
                       units               = '0.001'
                       scale_factor        = 0.001
                       add_offset          = 0
                       data_mode           = 'R'
                       ancillary_variables = 'PSAL_QC
I am a bit confused about what your code is supposed to do 
most of the data obtained with ncread are function of DEPTH andTIME and not longitude and latitude , so I am a bit sceptical about the fact that this is indeed suppose to give a value of for the Krigi lat / lon coordinates ! 
[~,ilatK]=min(abs(lat-3.4));% Indice latitude correspondant à Kribi 
[~,ilonK]=min(abs(lon-6.33));% Indice longitude correspondant à Kribi
if (~isempty(ilatK)&&~isempty(ilonK))
    sss=ncread([PathName,FileName],'PSAL',[ilonK ilatK],[1 1]);
    SSSK=NaN(1,size(sss,2)); % Initialisation du vecteur de sorties des pH 
    for i=1:size(sss,2)
        SSSK(i)=sss(:,:,i); % Compilation de la série temporelle
    end
    sal=ncread([PathName,FileName],'PSAL_QC',[ilonK ilatK],[1 1]);
    salK=NaN(1,size(sal,3)); % Initialisation du vecteur de sorties chl
    for i=1:size(sal,3)
        salK(i)=sal(:,:,i); % Compilation de la série temporelle
    end
end
I suppose that measurements are done by a device that measure parameters at one location at a time , so you don't have a map but a trajectory vs time . If you plot lat and lon vs time you get this : 

Weitere Antworten (1)
  Peter Perkins
    
 am 27 Nov. 2023
        Using datenums is no longer recomended.
t1=double(ncread([PathName,FileName],'TIME')); 
...
    t(i)=datenum(1900,1,1,0,0,0)+(t1(i)/24); % Conversion en jours juliens
...
    plot(t,pHK,'.-k','Linewidth',1,'Markersize',10) % Hauteur au large
If the time in the file is really the number of hours since 1900, stored as a number, do this:
t1 = ncread([PathName,FileName],'TIME'); 
...
    t(i) = datetime(1900,1,t1(i)); % Conversion en jours juliens
...
    plot(t,pHK,'.-k','Linewidth',1,'Markersize',10) % Hauteur au large
0 Kommentare
Siehe auch
Kategorien
				Mehr zu Data Import and Analysis finden Sie in Help Center und File Exchange
			
	Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!


