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

2 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

Categories

Find more on Programming 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!