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

3 views (last 30 days)
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')

Accepted Answer

Mathieu NOE
Mathieu NOE on 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 Comments
Mathieu NOE
Mathieu NOE on 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
Mathieu NOE on 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 :

Sign in to comment.

More Answers (1)

Peter Perkins
Peter Perkins on 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

Community Treasure Hunt

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

Start Hunting!