Gap in the plot
Hi,
I am given a contineuous water level data measured and measured cross section to calculate discharge by Manning equation. All these data were checked for gaps, and found contineuous. When running the following script the discharge having gaps, creats a jump which looks strange. Tried for several weeks to find it, but didn’t. Hope for your help.
%this script reads Cross and crocss-sections data from tributaries:
%Tzipori, Yiftachel, Tzvi, Oz, Mizra, Adashim, Taanach and Keini
%It 1. reads divers data 2. change them to be water level 3. calculate
%discharge.
%Discharge calculation requires:
%1. find the water level in the current time step 2. caculate discharge using manning
tic
clear all
close all
%%%%%%%%%%%% Part A: find water level of all sites%%%%%%%%%%%%%%%%%%%
% Get a list of all txt files in the current folder, or subfolders of it.
%Read discharge
lines1 = readlines("trib_names.txt");
trib_name=(categorical(lines1));
for i=1:length(trib_name)
trib1{i}=trib_name(i)
end
trib=string(trib1)+’level’; %for the case of water level;
cd ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischdivers’;
lvl=NaN(418907,16);%columns: longest diver data, rows/2= # of tributaries. each st. ha two columns: for date and data of stations
AllDivresfiles = natsortfiles(dir(‘*.txt’));
for k =1:size(AllDivresfiles,1)
% k=5%shut ‘end’ at the end of loop. Find thie in ‘kishon_catchment_discharge_calc_for1trib’
thisfilename1 = AllDivresfiles(k).name; %just the name
thisdata1 = load(thisfilename1); %load just this file
Lng_lvl(k)=length(thisdata1);
date1=x2mdate(thisdata1(:,1));
%plot(date1,thisdata1(:,2));ylabel(‘Level (cm)’)
num_days=length(unique(thisdata1(:,1)));
interval=720;%720 is the daily period of the time interval of data: 2 min. length(thisdata1(:,1))/(num_days*24);%Model has 1 hours intervals, so this is the interval obtained here. Obtain intervals according to 10 min intrval
baseP=min(thisdata1(:,2));%background pressure
curr_lvl1=(thisdata1(:,2)-baseP)/100;%units change to meter%for smoothing turn curr_lvl to curr_lvl1
curr_lvl=smoothdata(curr_lvl1,’movmean’,interval);
%curr_lvl=smoothdata(curr_lvl,’gaussian’,interval);
lvl(1:size(thisdata1,1),(2*k-1):2*k)=[date1 curr_lvl];%Removing diver’s bias (correct only forephemeral tributaries) and converting to m
lngt(k)=length(curr_lvl);
end
in0=find((lvl)==0);
lvl(in0)=NaN;
figure(‘Name’,’Adashim’);adash=lvl(:,1:2);plot(adash(:,1),adash(:,2));datetick;title=’Adashim lvl’;
figure(‘Name’,’Keini’);keini=lvl(:,3:4);plot(keini(:,1),keini(:,2));title=’Keini lvl’;datetick;
figure(‘Name’,’KishonMaale’);KishonMaale=lvl(:,5:6);plot(KishonMaale(:,1),KishonMaale(:,2));title=’KishonMaale lvl’;datetick
figure(‘Name’,’Mizra’);mizra=lvl(:,7:8);plot(mizra(:,1),mizra(:,2));title=’Mizra lvl’;datetick;
figure(‘Name’,’Oz’);oz=lvl(:,9:10);plot(oz(:,1),oz(:,2));title=’Oz lvl’;datetick;
figure(‘Name’,’Taanach’);taanach=lvl(:,11:12);plot(taanach(:,1),taanach(:,2));title=’Taanach lvl’;datetick;
figure(‘Name’,’Tzvi’);title=’Tzvi lvl’;tzvi=lvl(:,13:14);plot(tzvi(:,1),tzvi(:,2));datetick;
figure(‘Name’,’Yiftachel’);title=’Yiftachel lvl’;yiftachel=lvl(:,15:16);plot(yiftachel(:,1),yiftachel(:,2));title=’Yiftachel lvl’;datetick;
%%%%%%%%%%%% Part B: find cross section%%%%%%%%%%%%%%%%%%%
%Read cross sections. In most tributaries the divers are in the downstream
%All files with ‘1’ are upstrream and all with ‘2’ are downstreaצ. For example: OzUp=Oz1, OzDown=Oz2
trib=length(AllDivresfiles);
cd ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischcrossectionsDownstream’;
crs=NaN(max(lngt),5*trib);%columns: longest cross section data, rows-# of stations*5
AllCrossfiles = natsortfiles(dir(‘*.txt’));
min_for_slope1=NaN(trib,3);
%Collect all downstream crossection data together
for l = 1:length(AllCrossfiles)
%l=1 %shut ‘end’ at the end of loop
thisfilename2 = AllCrossfiles(l).name; %just the name
thisdata2 = load(thisfilename2); %load current file
crs(1:length(thisdata2),5*l-4:l*5)=thisdata2;
[MIN,I0]=min(thisdata2(:,3));%sea level elevation to have absolute elevation
min_for_slope1(l,:)=thisdata2(I0,1:3);%The slope is obtained by the slope bwtween two min of adjacent croo sections.
Lng(l)=length(thisdata2);
% figure(‘Name’,[char(trib_name(l)) ‘_crs_dn’])
% plot(crs(:,5*l-4),crs(:,5*l))
end
% adash1_crs=crs(:,4:5);keini1_crs=crs(:,9:10);mizra1_crs=crs(:,14:15);oz1_crs=crs(:,19:20);taanach1_crs=crs(:,24:25);tzvi1_crs=crs(:,29:30);yiftachel1_crs=crs(:,34:35);
% figure(11);plot(adash1_crs(:,1),adash1_crs(:,2));figure(22); plot(keini1_crs(:,1),keini1_crs(:,2));figure(33);plot(mizra1_crs(:,1),mizra1_crs(:,2));figure(44);plot(oz1_crs(:,1),oz1_crs(:,2));figure(55);plot(taanach1_crs(:,1),taanach1_crs(:,2));figure(66);plot(tzvi1_crs(:,1),tzvi1_crs(:,2));figure(77);plot(yiftachel1_crs(:,1),yiftachel1_crs(:,2));
%%%%% Part B1: Reading upstream for slope%%%%%%%%
min_for_slope2=NaN(trib,3);%coordinates (columns 1 and 2) and elevation (column3). of minimum in the cross sectioncolumns: longest diver data, rows/2= # of tributaries. each st. ha two columns: for date and data of stations
cd ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischcrossections’;
%min_for_slope=NaN(trib*2,2);%columns: longest cross section data, rows-# of stations
Allupstrmfiles = orderfields(dir(‘*.txt’));
for j = 1 : trib
upstrmName = Allupstrmfiles(j).name; %just the name; j or l
thisupstrm = load(upstrmName); %load just this file
[MIN1,I1]=min(thisupstrm(:,3));%Sea level to get absolute elevation
min_for_slope2(j,:)=(thisupstrm(I1,1:3));
% figure(‘Name’,[char(trib_name(j)) ‘_crs_up’])
% plot(thisupstrm(:,4),thisupstrm(:,3))
end
%%%%%%%%%%%%%%%%NEW 29052025%%%%%%%%%%%%%%%%%%%%%%%%%%%
all_lvl=NaN(max(lngt),trib);
discharge_all=NaN(418907,trib*2);%longest data
discharge=NaN(418907,2);
g=0;
f=0;
for t=1:trib
x=crs(1:Lng(t),5*t-4);
y=crs(1:Lng(t),5*t-3);
crsc=crs(1:Lng(t),5*t);
x0=x(1);y0=y(1);
distance1=((x-x0).^2+(y-y0).^2).^(1/2);%crs(1:Lng(t),5*t-1);
%Interpolation to create smoother discharge calculation
DF=abs(min(crsc(1:length(crsc)-1)-crsc(2:length(crsc))))/300;%Order of magnitude of distance difference and elveation are similar
basic_distance=distance1(1):DF:distance1(end);
distance11=1:(length(crsc));
distance=interp1(distance11,distance1,basic_distance);%to do interpolation to distance , since t is not a stright line
smth_crs=interp1(distance1,crsc,distance);
% smth_x=interp1(distance1,x,distance);
% smth_y=interp1(distance1,y,distance);
for m=1:Lng_lvl(t)%m=233648
if isnan(lvl(m,2*t))
discharge(m,1:2)=[lvl(m,2*t-1) NaN];
f=f+1;
continue
end
[M lvl_ind1]=min(abs(lvl(m,2*t)-smth_crs));%Find current water level in lvl and the current cross section
if smth_crs(lvl_ind1)==min(smth_crs)
discharge(m,1:2)=[lvl(m,2*t-1) 0];
g=g+1;
continue
end
lvl_ind2=find(smth_crs<=smth_crs(lvl_ind1));
distance2=distance(lvl_ind2);smth_crs2=smth_crs(lvl_ind2);
P1=((distance2(2:end)-distance2(1:end-1)).^2+(smth_crs2(2:end)-smth_crs2(1:end-1)).^2).^(1/2);
% lvl_ind2=find(smth_crs<=smth_crs(lvl_ind1));
% lvl_ind22=crsc<=min(crsc);
% crsc_lvl=crsc(lvl_ind22);
% x_lvl=x(lvl_ind22);%Due to wetted perimeter.
% y_lvl=y(lvl_ind22);%Due to wetted perimeter.
% % smth_x_lvl=smth_x(lvl_ind2);%Due to wetted perimeter.
% % smth_y_lvl=smth_y(lvl_ind2);%Due to wetted perimeter.
% smth_crs_lvl=smth_crs(lvl_ind2);%Due to wetted perimeter.
%Wetted perimeter:
% P1=((x_lvl(2:end)-x_lvl(1:end-1)).^2+(y_lvl(2:end)-y_lvl(1:end-1)).^2+(crsc_lvl(2:end)-crsc_lvl(1:end-1)).^2).^(1/2);
P(m)=sum(P1);
if isnan(P(m))
discharge(m,1:2)=[lvl(m,2*t-1) 0];
q=q+1;
continue
end
%Slope
dx = min_for_slope2(t,1) – min_for_slope1(t,1);
dy = min_for_slope2(t,2) – min_for_slope1(t,2);
dz(t) = min_for_slope2(t,3) – min_for_slope1(t,3);
distance3D = sqrt(dx^2 + dy^2 + dz(t)^2);
S(t) = abs(dz(t)) / distance3D;
% S(t)=min_for_slope2(t,3)-min_for_slope1(t,3)/(min_for_slope2(t,1)-min_for_slope1(t,1))^2+(min_for_slope2(t,2)-min_for_slope1(t,2))^2+…
% (min_for_slope2(t,3)-min_for_slope1(t,3))^2;%Slope. Constant for crossection
%Area
A(m) = trapz(distance(lvl_ind2), smth_crs(lvl_ind1) – smth_crs(lvl_ind2)); % assumes smth_crs is bed elevation
%A(m)=trapz(distance-smth_crs);
R(m)=A(m)/P(m);%Hydraulic radius;
%Manning
%Effect of water level on n of Manning:
if smth_crs(lvl_ind1) < 0.3
N = 0.1;
elseif and(smth_crs(lvl_ind1) >= 0.3, smth_crs(lvl_ind1) < 0.7)
N = 0.06;
else
N = 0.03;
end
discharge1=(A(m)*(1/N)*R(m)^(2/3))*S(t)^(1/2);
discharge(m,1:2)=[lvl(m,2*t-1) discharge1];
end%end ‘for m=…’
% figure(‘Name’,[char(trib_name(t)) ‘_smth_crs’])
% plot(smth_crs)
discharge_all(1:size(discharge,1),2*t-1:2*t)=discharge;
% [discharge2,ia,~] = unique(discharge1(:,2));
% discharge=[lvl(ia,2*t-1) discharge2];
% discharge_all(1:size(discharge,1),2*t-1:2*t)=discharge;
% plot(date1,discharge_all);
n=t;
discharge_all(1:length(discharge),2*n-1:2*n)=discharge;
% figure(n)
% plot(discharge_all(1:Lng_lvl(n),2*n-1),discharge_all(1:Lng_lvl(n),2*n))
% title=trib_name(n);
% ylabel(‘Q (m^3/s)’)
% datetick
figure(‘name’,char(trib_name(t)))
plot(lvl(:,2*t-1),discharge_all(:,2*t),’.’)
%title(names(t));
ylabel(‘Q (m^3/s)’)
datetick
end
%end %end for t=, m=…
discharge_all_smth=discharge_all;
filename = ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischdischarge_all_smth.mat’;
save( filename, ‘discharge_all_smth’ );
tocHi,
I am given a contineuous water level data measured and measured cross section to calculate discharge by Manning equation. All these data were checked for gaps, and found contineuous. When running the following script the discharge having gaps, creats a jump which looks strange. Tried for several weeks to find it, but didn’t. Hope for your help.
%this script reads Cross and crocss-sections data from tributaries:
%Tzipori, Yiftachel, Tzvi, Oz, Mizra, Adashim, Taanach and Keini
%It 1. reads divers data 2. change them to be water level 3. calculate
%discharge.
%Discharge calculation requires:
%1. find the water level in the current time step 2. caculate discharge using manning
tic
clear all
close all
%%%%%%%%%%%% Part A: find water level of all sites%%%%%%%%%%%%%%%%%%%
% Get a list of all txt files in the current folder, or subfolders of it.
%Read discharge
lines1 = readlines("trib_names.txt");
trib_name=(categorical(lines1));
for i=1:length(trib_name)
trib1{i}=trib_name(i)
end
trib=string(trib1)+’level’; %for the case of water level;
cd ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischdivers’;
lvl=NaN(418907,16);%columns: longest diver data, rows/2= # of tributaries. each st. ha two columns: for date and data of stations
AllDivresfiles = natsortfiles(dir(‘*.txt’));
for k =1:size(AllDivresfiles,1)
% k=5%shut ‘end’ at the end of loop. Find thie in ‘kishon_catchment_discharge_calc_for1trib’
thisfilename1 = AllDivresfiles(k).name; %just the name
thisdata1 = load(thisfilename1); %load just this file
Lng_lvl(k)=length(thisdata1);
date1=x2mdate(thisdata1(:,1));
%plot(date1,thisdata1(:,2));ylabel(‘Level (cm)’)
num_days=length(unique(thisdata1(:,1)));
interval=720;%720 is the daily period of the time interval of data: 2 min. length(thisdata1(:,1))/(num_days*24);%Model has 1 hours intervals, so this is the interval obtained here. Obtain intervals according to 10 min intrval
baseP=min(thisdata1(:,2));%background pressure
curr_lvl1=(thisdata1(:,2)-baseP)/100;%units change to meter%for smoothing turn curr_lvl to curr_lvl1
curr_lvl=smoothdata(curr_lvl1,’movmean’,interval);
%curr_lvl=smoothdata(curr_lvl,’gaussian’,interval);
lvl(1:size(thisdata1,1),(2*k-1):2*k)=[date1 curr_lvl];%Removing diver’s bias (correct only forephemeral tributaries) and converting to m
lngt(k)=length(curr_lvl);
end
in0=find((lvl)==0);
lvl(in0)=NaN;
figure(‘Name’,’Adashim’);adash=lvl(:,1:2);plot(adash(:,1),adash(:,2));datetick;title=’Adashim lvl’;
figure(‘Name’,’Keini’);keini=lvl(:,3:4);plot(keini(:,1),keini(:,2));title=’Keini lvl’;datetick;
figure(‘Name’,’KishonMaale’);KishonMaale=lvl(:,5:6);plot(KishonMaale(:,1),KishonMaale(:,2));title=’KishonMaale lvl’;datetick
figure(‘Name’,’Mizra’);mizra=lvl(:,7:8);plot(mizra(:,1),mizra(:,2));title=’Mizra lvl’;datetick;
figure(‘Name’,’Oz’);oz=lvl(:,9:10);plot(oz(:,1),oz(:,2));title=’Oz lvl’;datetick;
figure(‘Name’,’Taanach’);taanach=lvl(:,11:12);plot(taanach(:,1),taanach(:,2));title=’Taanach lvl’;datetick;
figure(‘Name’,’Tzvi’);title=’Tzvi lvl’;tzvi=lvl(:,13:14);plot(tzvi(:,1),tzvi(:,2));datetick;
figure(‘Name’,’Yiftachel’);title=’Yiftachel lvl’;yiftachel=lvl(:,15:16);plot(yiftachel(:,1),yiftachel(:,2));title=’Yiftachel lvl’;datetick;
%%%%%%%%%%%% Part B: find cross section%%%%%%%%%%%%%%%%%%%
%Read cross sections. In most tributaries the divers are in the downstream
%All files with ‘1’ are upstrream and all with ‘2’ are downstreaצ. For example: OzUp=Oz1, OzDown=Oz2
trib=length(AllDivresfiles);
cd ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischcrossectionsDownstream’;
crs=NaN(max(lngt),5*trib);%columns: longest cross section data, rows-# of stations*5
AllCrossfiles = natsortfiles(dir(‘*.txt’));
min_for_slope1=NaN(trib,3);
%Collect all downstream crossection data together
for l = 1:length(AllCrossfiles)
%l=1 %shut ‘end’ at the end of loop
thisfilename2 = AllCrossfiles(l).name; %just the name
thisdata2 = load(thisfilename2); %load current file
crs(1:length(thisdata2),5*l-4:l*5)=thisdata2;
[MIN,I0]=min(thisdata2(:,3));%sea level elevation to have absolute elevation
min_for_slope1(l,:)=thisdata2(I0,1:3);%The slope is obtained by the slope bwtween two min of adjacent croo sections.
Lng(l)=length(thisdata2);
% figure(‘Name’,[char(trib_name(l)) ‘_crs_dn’])
% plot(crs(:,5*l-4),crs(:,5*l))
end
% adash1_crs=crs(:,4:5);keini1_crs=crs(:,9:10);mizra1_crs=crs(:,14:15);oz1_crs=crs(:,19:20);taanach1_crs=crs(:,24:25);tzvi1_crs=crs(:,29:30);yiftachel1_crs=crs(:,34:35);
% figure(11);plot(adash1_crs(:,1),adash1_crs(:,2));figure(22); plot(keini1_crs(:,1),keini1_crs(:,2));figure(33);plot(mizra1_crs(:,1),mizra1_crs(:,2));figure(44);plot(oz1_crs(:,1),oz1_crs(:,2));figure(55);plot(taanach1_crs(:,1),taanach1_crs(:,2));figure(66);plot(tzvi1_crs(:,1),tzvi1_crs(:,2));figure(77);plot(yiftachel1_crs(:,1),yiftachel1_crs(:,2));
%%%%% Part B1: Reading upstream for slope%%%%%%%%
min_for_slope2=NaN(trib,3);%coordinates (columns 1 and 2) and elevation (column3). of minimum in the cross sectioncolumns: longest diver data, rows/2= # of tributaries. each st. ha two columns: for date and data of stations
cd ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischcrossections’;
%min_for_slope=NaN(trib*2,2);%columns: longest cross section data, rows-# of stations
Allupstrmfiles = orderfields(dir(‘*.txt’));
for j = 1 : trib
upstrmName = Allupstrmfiles(j).name; %just the name; j or l
thisupstrm = load(upstrmName); %load just this file
[MIN1,I1]=min(thisupstrm(:,3));%Sea level to get absolute elevation
min_for_slope2(j,:)=(thisupstrm(I1,1:3));
% figure(‘Name’,[char(trib_name(j)) ‘_crs_up’])
% plot(thisupstrm(:,4),thisupstrm(:,3))
end
%%%%%%%%%%%%%%%%NEW 29052025%%%%%%%%%%%%%%%%%%%%%%%%%%%
all_lvl=NaN(max(lngt),trib);
discharge_all=NaN(418907,trib*2);%longest data
discharge=NaN(418907,2);
g=0;
f=0;
for t=1:trib
x=crs(1:Lng(t),5*t-4);
y=crs(1:Lng(t),5*t-3);
crsc=crs(1:Lng(t),5*t);
x0=x(1);y0=y(1);
distance1=((x-x0).^2+(y-y0).^2).^(1/2);%crs(1:Lng(t),5*t-1);
%Interpolation to create smoother discharge calculation
DF=abs(min(crsc(1:length(crsc)-1)-crsc(2:length(crsc))))/300;%Order of magnitude of distance difference and elveation are similar
basic_distance=distance1(1):DF:distance1(end);
distance11=1:(length(crsc));
distance=interp1(distance11,distance1,basic_distance);%to do interpolation to distance , since t is not a stright line
smth_crs=interp1(distance1,crsc,distance);
% smth_x=interp1(distance1,x,distance);
% smth_y=interp1(distance1,y,distance);
for m=1:Lng_lvl(t)%m=233648
if isnan(lvl(m,2*t))
discharge(m,1:2)=[lvl(m,2*t-1) NaN];
f=f+1;
continue
end
[M lvl_ind1]=min(abs(lvl(m,2*t)-smth_crs));%Find current water level in lvl and the current cross section
if smth_crs(lvl_ind1)==min(smth_crs)
discharge(m,1:2)=[lvl(m,2*t-1) 0];
g=g+1;
continue
end
lvl_ind2=find(smth_crs<=smth_crs(lvl_ind1));
distance2=distance(lvl_ind2);smth_crs2=smth_crs(lvl_ind2);
P1=((distance2(2:end)-distance2(1:end-1)).^2+(smth_crs2(2:end)-smth_crs2(1:end-1)).^2).^(1/2);
% lvl_ind2=find(smth_crs<=smth_crs(lvl_ind1));
% lvl_ind22=crsc<=min(crsc);
% crsc_lvl=crsc(lvl_ind22);
% x_lvl=x(lvl_ind22);%Due to wetted perimeter.
% y_lvl=y(lvl_ind22);%Due to wetted perimeter.
% % smth_x_lvl=smth_x(lvl_ind2);%Due to wetted perimeter.
% % smth_y_lvl=smth_y(lvl_ind2);%Due to wetted perimeter.
% smth_crs_lvl=smth_crs(lvl_ind2);%Due to wetted perimeter.
%Wetted perimeter:
% P1=((x_lvl(2:end)-x_lvl(1:end-1)).^2+(y_lvl(2:end)-y_lvl(1:end-1)).^2+(crsc_lvl(2:end)-crsc_lvl(1:end-1)).^2).^(1/2);
P(m)=sum(P1);
if isnan(P(m))
discharge(m,1:2)=[lvl(m,2*t-1) 0];
q=q+1;
continue
end
%Slope
dx = min_for_slope2(t,1) – min_for_slope1(t,1);
dy = min_for_slope2(t,2) – min_for_slope1(t,2);
dz(t) = min_for_slope2(t,3) – min_for_slope1(t,3);
distance3D = sqrt(dx^2 + dy^2 + dz(t)^2);
S(t) = abs(dz(t)) / distance3D;
% S(t)=min_for_slope2(t,3)-min_for_slope1(t,3)/(min_for_slope2(t,1)-min_for_slope1(t,1))^2+(min_for_slope2(t,2)-min_for_slope1(t,2))^2+…
% (min_for_slope2(t,3)-min_for_slope1(t,3))^2;%Slope. Constant for crossection
%Area
A(m) = trapz(distance(lvl_ind2), smth_crs(lvl_ind1) – smth_crs(lvl_ind2)); % assumes smth_crs is bed elevation
%A(m)=trapz(distance-smth_crs);
R(m)=A(m)/P(m);%Hydraulic radius;
%Manning
%Effect of water level on n of Manning:
if smth_crs(lvl_ind1) < 0.3
N = 0.1;
elseif and(smth_crs(lvl_ind1) >= 0.3, smth_crs(lvl_ind1) < 0.7)
N = 0.06;
else
N = 0.03;
end
discharge1=(A(m)*(1/N)*R(m)^(2/3))*S(t)^(1/2);
discharge(m,1:2)=[lvl(m,2*t-1) discharge1];
end%end ‘for m=…’
% figure(‘Name’,[char(trib_name(t)) ‘_smth_crs’])
% plot(smth_crs)
discharge_all(1:size(discharge,1),2*t-1:2*t)=discharge;
% [discharge2,ia,~] = unique(discharge1(:,2));
% discharge=[lvl(ia,2*t-1) discharge2];
% discharge_all(1:size(discharge,1),2*t-1:2*t)=discharge;
% plot(date1,discharge_all);
n=t;
discharge_all(1:length(discharge),2*n-1:2*n)=discharge;
% figure(n)
% plot(discharge_all(1:Lng_lvl(n),2*n-1),discharge_all(1:Lng_lvl(n),2*n))
% title=trib_name(n);
% ylabel(‘Q (m^3/s)’)
% datetick
figure(‘name’,char(trib_name(t)))
plot(lvl(:,2*t-1),discharge_all(:,2*t),’.’)
%title(names(t));
ylabel(‘Q (m^3/s)’)
datetick
end
%end %end for t=, m=…
discharge_all_smth=discharge_all;
filename = ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischdischarge_all_smth.mat’;
save( filename, ‘discharge_all_smth’ );
toc Hi,
I am given a contineuous water level data measured and measured cross section to calculate discharge by Manning equation. All these data were checked for gaps, and found contineuous. When running the following script the discharge having gaps, creats a jump which looks strange. Tried for several weeks to find it, but didn’t. Hope for your help.
%this script reads Cross and crocss-sections data from tributaries:
%Tzipori, Yiftachel, Tzvi, Oz, Mizra, Adashim, Taanach and Keini
%It 1. reads divers data 2. change them to be water level 3. calculate
%discharge.
%Discharge calculation requires:
%1. find the water level in the current time step 2. caculate discharge using manning
tic
clear all
close all
%%%%%%%%%%%% Part A: find water level of all sites%%%%%%%%%%%%%%%%%%%
% Get a list of all txt files in the current folder, or subfolders of it.
%Read discharge
lines1 = readlines("trib_names.txt");
trib_name=(categorical(lines1));
for i=1:length(trib_name)
trib1{i}=trib_name(i)
end
trib=string(trib1)+’level’; %for the case of water level;
cd ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischdivers’;
lvl=NaN(418907,16);%columns: longest diver data, rows/2= # of tributaries. each st. ha two columns: for date and data of stations
AllDivresfiles = natsortfiles(dir(‘*.txt’));
for k =1:size(AllDivresfiles,1)
% k=5%shut ‘end’ at the end of loop. Find thie in ‘kishon_catchment_discharge_calc_for1trib’
thisfilename1 = AllDivresfiles(k).name; %just the name
thisdata1 = load(thisfilename1); %load just this file
Lng_lvl(k)=length(thisdata1);
date1=x2mdate(thisdata1(:,1));
%plot(date1,thisdata1(:,2));ylabel(‘Level (cm)’)
num_days=length(unique(thisdata1(:,1)));
interval=720;%720 is the daily period of the time interval of data: 2 min. length(thisdata1(:,1))/(num_days*24);%Model has 1 hours intervals, so this is the interval obtained here. Obtain intervals according to 10 min intrval
baseP=min(thisdata1(:,2));%background pressure
curr_lvl1=(thisdata1(:,2)-baseP)/100;%units change to meter%for smoothing turn curr_lvl to curr_lvl1
curr_lvl=smoothdata(curr_lvl1,’movmean’,interval);
%curr_lvl=smoothdata(curr_lvl,’gaussian’,interval);
lvl(1:size(thisdata1,1),(2*k-1):2*k)=[date1 curr_lvl];%Removing diver’s bias (correct only forephemeral tributaries) and converting to m
lngt(k)=length(curr_lvl);
end
in0=find((lvl)==0);
lvl(in0)=NaN;
figure(‘Name’,’Adashim’);adash=lvl(:,1:2);plot(adash(:,1),adash(:,2));datetick;title=’Adashim lvl’;
figure(‘Name’,’Keini’);keini=lvl(:,3:4);plot(keini(:,1),keini(:,2));title=’Keini lvl’;datetick;
figure(‘Name’,’KishonMaale’);KishonMaale=lvl(:,5:6);plot(KishonMaale(:,1),KishonMaale(:,2));title=’KishonMaale lvl’;datetick
figure(‘Name’,’Mizra’);mizra=lvl(:,7:8);plot(mizra(:,1),mizra(:,2));title=’Mizra lvl’;datetick;
figure(‘Name’,’Oz’);oz=lvl(:,9:10);plot(oz(:,1),oz(:,2));title=’Oz lvl’;datetick;
figure(‘Name’,’Taanach’);taanach=lvl(:,11:12);plot(taanach(:,1),taanach(:,2));title=’Taanach lvl’;datetick;
figure(‘Name’,’Tzvi’);title=’Tzvi lvl’;tzvi=lvl(:,13:14);plot(tzvi(:,1),tzvi(:,2));datetick;
figure(‘Name’,’Yiftachel’);title=’Yiftachel lvl’;yiftachel=lvl(:,15:16);plot(yiftachel(:,1),yiftachel(:,2));title=’Yiftachel lvl’;datetick;
%%%%%%%%%%%% Part B: find cross section%%%%%%%%%%%%%%%%%%%
%Read cross sections. In most tributaries the divers are in the downstream
%All files with ‘1’ are upstrream and all with ‘2’ are downstreaצ. For example: OzUp=Oz1, OzDown=Oz2
trib=length(AllDivresfiles);
cd ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischcrossectionsDownstream’;
crs=NaN(max(lngt),5*trib);%columns: longest cross section data, rows-# of stations*5
AllCrossfiles = natsortfiles(dir(‘*.txt’));
min_for_slope1=NaN(trib,3);
%Collect all downstream crossection data together
for l = 1:length(AllCrossfiles)
%l=1 %shut ‘end’ at the end of loop
thisfilename2 = AllCrossfiles(l).name; %just the name
thisdata2 = load(thisfilename2); %load current file
crs(1:length(thisdata2),5*l-4:l*5)=thisdata2;
[MIN,I0]=min(thisdata2(:,3));%sea level elevation to have absolute elevation
min_for_slope1(l,:)=thisdata2(I0,1:3);%The slope is obtained by the slope bwtween two min of adjacent croo sections.
Lng(l)=length(thisdata2);
% figure(‘Name’,[char(trib_name(l)) ‘_crs_dn’])
% plot(crs(:,5*l-4),crs(:,5*l))
end
% adash1_crs=crs(:,4:5);keini1_crs=crs(:,9:10);mizra1_crs=crs(:,14:15);oz1_crs=crs(:,19:20);taanach1_crs=crs(:,24:25);tzvi1_crs=crs(:,29:30);yiftachel1_crs=crs(:,34:35);
% figure(11);plot(adash1_crs(:,1),adash1_crs(:,2));figure(22); plot(keini1_crs(:,1),keini1_crs(:,2));figure(33);plot(mizra1_crs(:,1),mizra1_crs(:,2));figure(44);plot(oz1_crs(:,1),oz1_crs(:,2));figure(55);plot(taanach1_crs(:,1),taanach1_crs(:,2));figure(66);plot(tzvi1_crs(:,1),tzvi1_crs(:,2));figure(77);plot(yiftachel1_crs(:,1),yiftachel1_crs(:,2));
%%%%% Part B1: Reading upstream for slope%%%%%%%%
min_for_slope2=NaN(trib,3);%coordinates (columns 1 and 2) and elevation (column3). of minimum in the cross sectioncolumns: longest diver data, rows/2= # of tributaries. each st. ha two columns: for date and data of stations
cd ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischcrossections’;
%min_for_slope=NaN(trib*2,2);%columns: longest cross section data, rows-# of stations
Allupstrmfiles = orderfields(dir(‘*.txt’));
for j = 1 : trib
upstrmName = Allupstrmfiles(j).name; %just the name; j or l
thisupstrm = load(upstrmName); %load just this file
[MIN1,I1]=min(thisupstrm(:,3));%Sea level to get absolute elevation
min_for_slope2(j,:)=(thisupstrm(I1,1:3));
% figure(‘Name’,[char(trib_name(j)) ‘_crs_up’])
% plot(thisupstrm(:,4),thisupstrm(:,3))
end
%%%%%%%%%%%%%%%%NEW 29052025%%%%%%%%%%%%%%%%%%%%%%%%%%%
all_lvl=NaN(max(lngt),trib);
discharge_all=NaN(418907,trib*2);%longest data
discharge=NaN(418907,2);
g=0;
f=0;
for t=1:trib
x=crs(1:Lng(t),5*t-4);
y=crs(1:Lng(t),5*t-3);
crsc=crs(1:Lng(t),5*t);
x0=x(1);y0=y(1);
distance1=((x-x0).^2+(y-y0).^2).^(1/2);%crs(1:Lng(t),5*t-1);
%Interpolation to create smoother discharge calculation
DF=abs(min(crsc(1:length(crsc)-1)-crsc(2:length(crsc))))/300;%Order of magnitude of distance difference and elveation are similar
basic_distance=distance1(1):DF:distance1(end);
distance11=1:(length(crsc));
distance=interp1(distance11,distance1,basic_distance);%to do interpolation to distance , since t is not a stright line
smth_crs=interp1(distance1,crsc,distance);
% smth_x=interp1(distance1,x,distance);
% smth_y=interp1(distance1,y,distance);
for m=1:Lng_lvl(t)%m=233648
if isnan(lvl(m,2*t))
discharge(m,1:2)=[lvl(m,2*t-1) NaN];
f=f+1;
continue
end
[M lvl_ind1]=min(abs(lvl(m,2*t)-smth_crs));%Find current water level in lvl and the current cross section
if smth_crs(lvl_ind1)==min(smth_crs)
discharge(m,1:2)=[lvl(m,2*t-1) 0];
g=g+1;
continue
end
lvl_ind2=find(smth_crs<=smth_crs(lvl_ind1));
distance2=distance(lvl_ind2);smth_crs2=smth_crs(lvl_ind2);
P1=((distance2(2:end)-distance2(1:end-1)).^2+(smth_crs2(2:end)-smth_crs2(1:end-1)).^2).^(1/2);
% lvl_ind2=find(smth_crs<=smth_crs(lvl_ind1));
% lvl_ind22=crsc<=min(crsc);
% crsc_lvl=crsc(lvl_ind22);
% x_lvl=x(lvl_ind22);%Due to wetted perimeter.
% y_lvl=y(lvl_ind22);%Due to wetted perimeter.
% % smth_x_lvl=smth_x(lvl_ind2);%Due to wetted perimeter.
% % smth_y_lvl=smth_y(lvl_ind2);%Due to wetted perimeter.
% smth_crs_lvl=smth_crs(lvl_ind2);%Due to wetted perimeter.
%Wetted perimeter:
% P1=((x_lvl(2:end)-x_lvl(1:end-1)).^2+(y_lvl(2:end)-y_lvl(1:end-1)).^2+(crsc_lvl(2:end)-crsc_lvl(1:end-1)).^2).^(1/2);
P(m)=sum(P1);
if isnan(P(m))
discharge(m,1:2)=[lvl(m,2*t-1) 0];
q=q+1;
continue
end
%Slope
dx = min_for_slope2(t,1) – min_for_slope1(t,1);
dy = min_for_slope2(t,2) – min_for_slope1(t,2);
dz(t) = min_for_slope2(t,3) – min_for_slope1(t,3);
distance3D = sqrt(dx^2 + dy^2 + dz(t)^2);
S(t) = abs(dz(t)) / distance3D;
% S(t)=min_for_slope2(t,3)-min_for_slope1(t,3)/(min_for_slope2(t,1)-min_for_slope1(t,1))^2+(min_for_slope2(t,2)-min_for_slope1(t,2))^2+…
% (min_for_slope2(t,3)-min_for_slope1(t,3))^2;%Slope. Constant for crossection
%Area
A(m) = trapz(distance(lvl_ind2), smth_crs(lvl_ind1) – smth_crs(lvl_ind2)); % assumes smth_crs is bed elevation
%A(m)=trapz(distance-smth_crs);
R(m)=A(m)/P(m);%Hydraulic radius;
%Manning
%Effect of water level on n of Manning:
if smth_crs(lvl_ind1) < 0.3
N = 0.1;
elseif and(smth_crs(lvl_ind1) >= 0.3, smth_crs(lvl_ind1) < 0.7)
N = 0.06;
else
N = 0.03;
end
discharge1=(A(m)*(1/N)*R(m)^(2/3))*S(t)^(1/2);
discharge(m,1:2)=[lvl(m,2*t-1) discharge1];
end%end ‘for m=…’
% figure(‘Name’,[char(trib_name(t)) ‘_smth_crs’])
% plot(smth_crs)
discharge_all(1:size(discharge,1),2*t-1:2*t)=discharge;
% [discharge2,ia,~] = unique(discharge1(:,2));
% discharge=[lvl(ia,2*t-1) discharge2];
% discharge_all(1:size(discharge,1),2*t-1:2*t)=discharge;
% plot(date1,discharge_all);
n=t;
discharge_all(1:length(discharge),2*n-1:2*n)=discharge;
% figure(n)
% plot(discharge_all(1:Lng_lvl(n),2*n-1),discharge_all(1:Lng_lvl(n),2*n))
% title=trib_name(n);
% ylabel(‘Q (m^3/s)’)
% datetick
figure(‘name’,char(trib_name(t)))
plot(lvl(:,2*t-1),discharge_all(:,2*t),’.’)
%title(names(t));
ylabel(‘Q (m^3/s)’)
datetick
end
%end %end for t=, m=…
discharge_all_smth=discharge_all;
filename = ‘C:UsersuserDocumentsShulamitPhDMatlab_PhDHydrologyDischdischarge_all_smth.mat’;
save( filename, ‘discharge_all_smth’ );
toc hydrology manning water level gap in discharge MATLAB Answers — New Questions
