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");

Error using readlines (line 72)
Unable to find or open 'trib_names.txt'. Check the path and filename or file permissions.

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:\Users\user\Documents\Shulamit\PhD\Matlab_PhD\Hydrology\Disch\divers';
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:\Users\user\Documents\Shulamit\PhD\Matlab_PhD\Hydrology\Disch\crossections\Downstream';
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:\Users\user\Documents\Shulamit\PhD\Matlab_PhD\Hydrology\Disch\crossections';

%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:\Users\user\Documents\Shulamit\PhD\Matlab_PhD\Hydrology\Disch\discharge_all_smth.mat';
save( filename, 'discharge_all_smth' );
toc

[Steven Lord]

We don't have your data so we can't run this code and see what the picture looks like. We also can't see the picture generated on your machine. Please either attach a small data set that reproduces the problem or a .fig file / image file demonstrating the problem.

If I had to guess, I'd say that calling allfinite on the data you're plotting is going to return false. The plot function doesn't display a point if either of its coordinates is a missing value like NaN. If it does, you can use ismissing and find to identify the missing values.

x = [1:5 NaN 7:10];
areAnyValuesMissing = allfinite(x)

areAnyValuesMissing = logical
   0

missingMask = ismissing(x)

missingMask = 1×10 logical array
   0   0   0   0   0   1   0   0   0   0

missingLocations = find(missingMask)

missingLocations = 6

plot(x, 1:10, '-o') % Note that the line is broken at x = 6
grid on

●

[Stephen23]

Quite possibly the plotted data have NaN values in them: given that you call the NaN function ten times, NaN values seem very likely.

Also: avoid CD in code. Using absolute/relative filenames is more efficient and makes debugging much easier.

[Sam Chak]

Hi @Shulamit Nussboim, The "trib_names.txt" file is not provided. Consequently, other users in the forum are unable to run the code and investigate the issue on this platform.

Answer 1

[Shulamit Nussboim]

OK, thank you all. The problem was solved. It was since the options of n of Manning (if loop:0.01,0.03 and 0.06). The jump from intermediate to extremely low n resulted in sudden lowering of discharge, resulted in the jump.

(The NaNs you proposed to check is not profitted the jump location in data).