-
Notifications
You must be signed in to change notification settings - Fork 35
/
JunctionAngle.m
1497 lines (1331 loc) · 44.3 KB
/
JunctionAngle.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function [junctions,IX,varargout]=JunctionAngle(S,A,DEM,fit_distance,varargin);
% Usage:
% [junctions,IX]=JunctionAngle(S,A,DEM,fit_distance);
% [junctions,IX]=JunctionAngle(S,A,DEM,fit_distance,'name',value);
% [junctions,IX,mean_junctions]=JunctionAngle(S,A,DEM,[# # #]);
%
% Description:
% Function calculates the angle of stream junctions. The basic strategy and nomenclature
% used in this function is adapted from Howard, 1971, 'Optimal angles of stream junction:
% Geometric, stability to capture, and minimum power criteria'. The code first fits the
% upstream and downstream links that meet at junctions with linear fits. The lengths of
% stream used to perform these fits is controlled with the 'fit_distance' parameter. Junction
% angles are calculated as a pair of angles (e1 and e2) between each tributary (tributary 1 and
% tributary 2) and the upstream projection of the downstream link. The total junction angle
% is calculated as the angle between the two upstream links, which in well behaved junctions
% should be the sum of e1 and e2, but is not always the case (i.e. for some junctions both
% upstream links are on the same side of the upstream projection of the downstream link, in these
% cases the junction angle will not be the sum of the e1 and e2 angles, but instead will be the true
% angle between link 1 and link 2, whereas the values of e1 ane e2 will be measured from the upstream plane).
% The function also calculates predicted junction angles based on the simple geometric criteria
% described in Howard, 1971, and the handedness of junction angle as defined by James & Krumbein,
% 1969, 'Frequency distributions of stream link lengths'. It is important to note that angles
% related to junctions at which more than two upstream links meet or for which there are no
% downstream nodes (i.e. the junction is an outlet) are set to NaN as these are undefined in
% the Howard criteria. If a GRIDobj of estimated discharges in m^3/sec is provided, then an
% alternate prediction method using the minimum power criteria described in Howard, 1971 will
% also be calculated.
%
% Required Inputs:
% S - STREAMobj for calculating junction angles
% A - GRIDobj of flow accumulation grid
% DEM - GRIDobj of elevations, used for calculating slopes
% fit_distance - distance in stream distance (map units) for fitting stream links,
% if more than one value is provided to fit_distance, it is assumed you want to calculate
% junction angles fitting stream segments with the range of distances. The provided value(s)
% can be interpeted as the number of stream nodes over which to fit the stream links if the
% optional 'use_n_nodes' is set to true
%
% Optional Inputs:
% use_n_nodes [false] - logical flag to interpret the value(s) provided to fit_distance as
% the number of stream nodes to extract up and downstream as opposed to streamwise distance
% previous_IX [] - Calculation of the 'IX' output is the most time consuming aspect of the
% calculation, so if you wish to rerun the function with different parameters (e.g.
% recalculating predicted angles with a different reference concavity or different method)
% the presupplying the IX result from a previous run can be useful. It is important that
% this IX result was produced in a run using the same S, A, and DEM inputs and that the maximum
% value provided to fit distance does not exceed the maximum fit distance that was provided when
% you initially ran the function to produce the provided IX input. The code will check that the fit
% distances are compatible, but will not explicity check that S, A, and the DEM are the same.
% discharge [] - GRIDobj of discharge in m^3/sec to use in the area or stream power junction angle
% estimates. The most direct method for providing this would be using the FLOWobj and a GRIDobj
% of mean precipitation in m/sec in the 'flowacc' function to produce an estimate of discharge.
% If an argument is provided to this parameter, then this value will be used in place of area in
% the 'area' prediction method. This argument is required if the 'predict_angle_method' is
% 'minimum_power'.
% predict_angle_method ['area'] - method to use to calculate predicted junction angles. Both
% methods are simple geometric predictions from Howard, 1971. Valid inputs are 'slope',
% 'area', 'minimum_power', or 'all'. If pred_angle_method is set to 'area' (default), then you should
% make sure the ref_concavity is valid for the network you are analyzing, default value is 0.5. For the 'slope'
% method, the slope of the two upstream and one downstream link will be the mean of the gradient
% of those links over the length of stream sampled (controlled by the values provided to fit_distance).
% For the 'minimum_power' method, an argument is required for 'discharge'. If an entry for 'discharge' is provided
% then the 'area' method will also use this estimated discharge in place of a drainage area approximation.
% Providing 'all' will calculate both the slope and area predicted angles if no discharge is provided and
% will calculate the slope, area, and minimum power angles if a discharge is provided.
% ref_concavity - [0.5] - reference concavity for calculating predicted junction angles using
% Howard, 1971 area relation. If no value is provided, default value of 0.5 will be used.
% make_shape [false] - logical flag to generate a shapefile containing the junction locations
% and information stored in the various outputs
% verbose [false] - logical flag to report progress through function. This can be useful because
% certain steps of the process are very time consuming.
%
% Outputs:
% junctions - table (if single entry is provided to fit_distance) or a cell array
% where the first row contains the fit distances and the second row contains the
% table output for that fit distance. Columns of the table are:
%
% junction_number - ID number of junction
% junction_x - x coordinate of junction
% junction_y - y coordinate of junction
% junction_angle - angle in degrees between tributary 1 (e1) and tributary 2 (e2)
% handedness - handedness of junction as defined by James & Krumbein, 1969, options
% are Right, Left, or Undefined, stored as a categorical
% split - logical indicating whether upstream projection of downstream link lies between
% tributary 1 and tributary 2 (true) or does not (false)
% e1_obs_angle - angle between tributary 1 and upstream projection of downstream link
% e1_Apred_angle - predicted angle between tributary 1 and upstream projection of downstream link
% based on Howard 1971 area relations (will appear if 'predict_angle_method' is 'area' or 'all')
% e1_Spred_angle - predicted angle between tributary 1 and upstream projection of downstream link
% based on Howard 1971 slope relations (will appear if 'predict_angle_method' is 'slope' or 'all')
% e1_MPpred_angle - predicted angle between tributary 1 and upstream projection of downstream link
% based on Howard 1971 minimum power criteria (will appear if 'predict_angle_method' is 'minimum_power' or 'all'),
% but requires that a valid GRIDobj is provided to 'discharge'
% e1_rotation - orientation of tributary 1 with respect to the upstream projection of the
% downstream link, options are CCW (counter-clockwise), CW (clockwise), or Undefined
% (more than two tributaries meet at junction). If split is true, one of the upstream links will be
% CW and one will be CCW.
% e1_shreve - shreve order of tributary 1
% e1_direction - flow azimuth of tributary 1
% e1_distance - stream distance used to fit tributary 1
% e1_num - number of upstream nodes used to fit tributary 1
% e1_R2 - R2 value on linear fit of tributary 1
% e2_obs_angle - angle between tributary 2 and upstream projection of downstream link
% e2_Apred_angle - predicted angle between tributary 2 and upstream projection of downstream link
% based on Howard 1971 area relations (will appear if 'predict_angle_method' is 'area' or 'all')
% e2_Spred_angle - predicted angle between tributary 2 and upstream projection of downstream link
% based on Howard 1971 slope relations (will appear if 'predict_angle_method' is 'slope' or 'all')
% e2_MPpred_angle - predicted angle between tributary 1 and upstream projection of downstream link
% based on Howard 1971 minimum power criteria (will appear if 'predict_angle_method' is 'minimum_power' or 'all'),
% but requires that a valid GRIDobj is provided to 'discharge'
% e2_rotation - orientation of tributary 2 with respect to the upstream projection of the
% downstream link, options are CCW (counter-clockwise), CW (clockwise), or Undefined
% (more than two tributaries meet at junction)
% e2_shreve - shreve order of tributary 2
% e2_direction - flow azimuth of tributary 2
% e2_distance - stream distance used to fit tributary 2
% e2_num - number of upstream nodes used to fit tributary 2
% e2_R2 - R2 value on linear fit of tributary 2
% eS_direction - flow azimuth of downstream link
% eS_distance - stream distance used to fit downstream link
% eS_num - number of downstream nodes used to fit downstream link
% eS_R2 - R2 value of linear fit of downstream link
% IX - nx4 cell array with columns containing indices for selected positions
% in the node attributed list of the provided STREAMobj. Columns are confluence
% point, the selected downstream links, selected upstream links for tributary 1,
% and selected upstream links for tributary 2. If multiple values are provided to fit_distance
% the indices contained in the IX output will be those for the maximum fit_distance. Indices for
% nodes for upstream and downstream nodes are in topological order, i.e. the first index is the first
% upstream or downstream node and so on.
% mean_junctions - if multiple values are provided for fit distance, there is a third optional output
% which will be a table containing mean and standard deviations across the angles calculated using
% the different fit distances, including the mean total junction angle, standard deviation of
% the total junction angle, mean angle between tributary 1 and the upstream projection of the downstream
% link (e1), standard deviation of e1, mean angle between tributary 2 and the upstream projection of the
% downstream link (e2), standard deviation of e2. If the predict_angle_method is set to either 'slope'
% or 'both', this will also include the mean and standard deviations of the predicted e1 and e2 angles
% based on the slope method. Means and standard deviations are not calculated for the area method because
% these will not vary as a function of fit distance.
%
% Note:
% For all outputs, the tributaries are sorted so that tributary 1 is always the 'larger' of the
% two tributaries. The size of the tributaries are based on the shreve order of the two tributaries.
% In the event that the two tributaries have the same shreve order then the drainage area is used to
% determine respective size.
%
% Examples:
% [J,IX]=JunctionAngle(S,A,DEM,1000);
% [J,IX]=JunctionAngle(S,A,DEM,1000,'verbose',true,'make_shape',true,'ref_concavity',0.45);
% [J,IX,MJ]=JunctionAngle(S,A,[250 500 1000 2500 5000]);
% [J,IX]=JunctionAngle(S,A,DEM,1000,'predict_angle_method','slope');
% % Example for minimum power criteria
% PRECIPr=resample(PRECIP_M_S,DEM,'nearest');
% Q=flowacc(FD,PRECIPr);
% [J,IX]=JunctionAngle(S,A,DEM,1000,'discharge',Q,'predict_angle_method','minimum_power');
%
% Related Functions:
% InspectJunction JunctionLinks JunctionSegments
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function Written by Adam M. Forte - Updated : 09/27/21 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parse Inputs
p = inputParser;
p.FunctionName = 'JunctionAngle';
addRequired(p,'S',@(x) isa(x,'STREAMobj'));
addRequired(p,'A',@(x) isa(x,'GRIDobj'));
addRequired(p,'DEM',@(x) isa(x,'GRIDobj'));
addRequired(p,'fit_distance',@(x) isnumeric(x));
addParameter(p,'use_n_nodes',false,@(x) isscalar(x) && islogical(x));
addParameter(p,'previous_IX',[],@(x) iscell(x) || isempty(x));
addParameter(p,'discharge',[],@(x) isa(x,'GRIDobj'));
addParameter(p,'predict_angle_method','area',@(x) ischar(validatestring(x,{'slope','area','minimum_power','all'})));
addParameter(p,'ref_concavity',0.50,@(x) isscalar(x) && isnumeric(x));
addParameter(p,'make_shape',false,@(x) isscalar(x) && islogical(x));
addParameter(p,'verbose',false,@(x) isscalar(x) && islogical(x));
addParameter(p,'out_dir',[],@(x) isdir(x));
addParameter(p,'out_name',[],@(x) ischar(x));
parse(p,S,A,DEM,fit_distance,varargin{:});
S=p.Results.S;
A=p.Results.A;
DEM=p.Results.DEM;
fit_distance=p.Results.fit_distance;
use_n_nodes=p.Results.use_n_nodes;
PIX=p.Results.previous_IX;
Q=p.Results.discharge;
method=p.Results.predict_angle_method;
ref_theta=p.Results.ref_concavity;
make_shape=p.Results.make_shape;
verbose=p.Results.verbose;
out_dir=p.Results.out_dir;
out_name=p.Results.out_name;
if isempty(out_dir)
out_dir=pwd;
end
if use_n_nodes
if numel(fit_distance)==1
multi=false;
n=fit_distance;
if n<1
n=1;
end
fit_distance=n*(hypot(S.cellsize,S.cellsize));
else
multi=true;
n=max(fit_distance);
if n<1
n=1;
end
nlist=fit_distance;
nlist(nlist<1)=1;
nlist=unique(nlist);
if numel(nlist)==1
n=nlist;
multi=false;
end
fit_distance=nlist.*(hypot(S.cellsize,S.cellsize));
fit_distance=sort(fit_distance);
end
else
% Convert distance to number of nodes and check if multiple distances were provided.
if numel(fit_distance)==1
multi=false;
n=floor(fit_distance/hypot(S.cellsize,S.cellsize));
if n<1
if isdeployed
warndlg('Input fit_distance is too short, will not average across more than one node, setting to minimum distance')
end
warning('Input fit_distance is too short, will not average across more than one node, setting to minimum distance');
n=1;
end
else
n=floor(max(fit_distance)/hypot(S.cellsize,S.cellsize));
multi=true;
if n<1
if isdeployed
warndlg('All input fit_distance are too short, will not average across more than one node, setting to single minimum distance')
end
warning('All input fit_distance are too short, will not average across more than one node, setting to single minimum distance');
n=1;
mutli=false;
end
nlist=floor(fit_distance/hypot(S.cellsize,S.cellsize));
if any(nlist)<1
if isdeployed
warndlg('Some input fit_distance are too short, setting these to minimum distance')
end
warning('Some input fit_distance are too short, setting these to minimum distance');
nlist(nlist<1)=1;
nlist=unique(nlist);
if numel(nlist)==1
n=nlist;
multi=false;
end
end
fit_distance=sort(fit_distance);
end
end
% Extract stream distances
dst=S.distance;
% Calculate Shreve stream orders
so=streamorder(S,'shreve');
% Check for required inputs
if strcmp(method,'minimum_power') & isempty(Q)
error('Must provide a GRIDobj for the discharge parameter to use the minimum power method');
end
if ~isempty(Q) & ~validatealignment(S,Q)
error('STREAMobj and GRIDobj provided to discharge do not appear to have been derived from the same FLOWobj');
end
% Calculate gradient if pred_angle_method is slope
switch method
case {'slope','minimum_power','all'}
if verbose
disp('Calculating stream gradients')
end
z=mincosthydrocon(S,DEM,'interp',0.1);
sg=gradient(S,z);
G=GRIDobj(DEM);
G.Z(S.IXgrid)=sg;
end
% Find indices of nodes up and downstream
if isempty(PIX)
IX=nnodesupdown(S,A,n,so,verbose);
else
% Check that max distance is consistent
UPIX=PIX(:,2);
nn=cellfun(@numel,UPIX);
maxn=max(nn);
if maxn<n
if isdeployed
errordlg('The provided input to "previous_IX" is incompatible with the maximum value within "fit_distance"')
end
error('The provided input to "previous_IX" is incompatible with the maximum value within "fit_distance"');
else
IX=PIX;
end
end
% Convert stream x,y coordinates to lat-lon to preserve angular relationships, unless
% there is no projection (e.g. for model results) or projection is unrecognized
if verbose
disp('Converting to geographic coordinates...')
end
if isempty(S.georef)
xl=S.x; yl=S.y;
if isdeployed
warndlg('No projection was found, angles will be calculated based on projected coordinates')
end
warning('No projection was found, angles will be calculated based on projected coordinates')
projcoord=false;
else
try
[yl,xl]=projinv(S.georef.mstruct,S.x,S.y);
projcoord=true;
catch
xl=S.x; yl=S.y;
if isdeployed
warndlg('Projection was not recognized, angles will be calculated based on projected coordinates')
end
warning('Projection was not recognized, angles will be calculated based on projected coordinates')
projcoord=false;
end
end
% Break up output of node indicees
num_con=size(IX,1);
con=IX(:,1);
ds=IX(:,2);
us=IX(:,3:end);
% Check and mark if there any streams with more than 2 tribs and if any
% downstream segments of confluences don't exist to populate logical index
if size(us,2)>2
usvld=cellfun(@isempty,us(:,3));
else
usvld=true(size(con));
end
dsvld=~cellfun(@isempty,ds);
vld=dsvld & usvld;
if verbose
disp('Finding junction angles...')
end
if ~multi
% Preallocate arrays
link_angle=zeros(num_con,11);
strm_dir=zeros(num_con,3);
fit_streams=zeros(num_con,9);
switch method
case {'area','slope'}
pred_angle=zeros(num_con,2);
case 'both'
pred_angle=zeros(num_con,4);
end
for ii=1:num_con
if vld(ii)
% Extract x y node lists
% Confluence point
xc=xl(con{ii}); yc=yl(con{ii});
% Upstream links
xu1=xl(us{ii,1}); yu1=yl(us{ii,1});
xu2=xl(us{ii,2}); yu2=yl(us{ii,2});
% Downstream link
xd=xl(ds{ii,1}); yd=yl(ds{ii,1});
% Translate so all links originate from confluence
xu1=xu1-xc; xu2=xu2-xc;
yu1=yu1-yc; yu2=yu2-yc;
xd=xd-xc; yd=yd-yc;
%% Depending on orientation of links, doing a linear fit
% on the link may produce spurious results. To improve
% the result, each link is rotated towards horizontal,
% using the mean orientation of the link to find the angle
% of rotation
% Find mean angle of stream points from stream orientation
theta1ap=median(atan2(yu1,xu1));
theta2ap=median(atan2(yu2,xu2));
theta3ap=median(atan2(yd,xd));
% Extract and calculate stream distances
e1_dst=max(dst(us{ii,1}))-dst(con{ii});
e2_dst=max(dst(us{ii,2}))-dst(con{ii});
ds_dst=dst(con{ii})-min(dst(ds{ii,1}));
% Calculate predicted angles
switch method
case 'area'
% Howard 1971 area method
ix1=S.IXgrid(us{ii,1});
ix2=S.IXgrid(us{ii,2});
if isempty(Q)
da1=max(A.Z(ix1).*A.cellsize^2);
da2=max(A.Z(ix2).*A.cellsize^2);
else
da1=max(Q.Z(ix1));
da2=max(Q.Z(ix2));
end
e1p=acos(((da1+da2)/da1)^-ref_theta);
e2p=acos(((da1+da2)/da2)^-ref_theta);
% Convert to degrees
e1p=rad2deg(e1p);
e2p=rad2deg(e2p);
% Store out
pred_angle(ii,:)=[e1p e2p];
case 'slope'
% Howard 1971 slope method
ix1=S.IXgrid(us{ii,1});
ix2=S.IXgrid(us{ii,2});
ix3=S.IXgrid(ds{ii,1});
sl1=mean(G.Z(ix1));
sl2=mean(G.Z(ix2));
sl3=mean(G.Z(ix3));
r1=sl3/sl1;
r2=sl3/sl2;
if r1>1
e1p=0;
else
e1p=acos(r1);
end
if r2>1
e2p=0;
else
e2p=acos(r2);
end
% Convert to degrees
e1p=rad2deg(e1p);
e2p=rad2deg(e2p);
% Store out
pred_angle(ii,:)=[e1p e2p];
case 'minimum_power'
% Howard 1971 minimum power method
ix1=S.IXgrid(us{ii,1});
ix2=S.IXgrid(us{ii,2});
ix3=S.IXgrid(ds{ii,1});
sl1=mean(G.Z(ix1));
sl2=mean(G.Z(ix2));
sl3=mean(G.Z(ix3));
q1=max(Q.Z(ix1));
q2=max(Q.Z(ix2));
q3=q1+q2;
C1=q1*sl1;
C2=q2*sl2;
C3=q3*sl3;
inner1=((C3^2)+(C1^2)-(C2^2))/(2*C1*C3);
inner2=((C3^2)+(C2^2)-(C1^2))/(2*C2*C3);
if inner1>1
e1p=0;
else
e1p=acos(inner1);
end
if inner>1
e2p=0;
else
e2p=acos(inner2);
end
% Convert to degrees
e1p=rad2deg(e1p);
e2p=rad2deg(e2p);
% Store out
pred_angle(ii,:)=[e1p e2p];
case 'all'
% Howard 1971 area method
ix1=S.IXgrid(us{ii,1});
ix2=S.IXgrid(us{ii,2});
if isempty(Q)
da1=max(A.Z(ix1).*A.cellsize^2);
da2=max(A.Z(ix2).*A.cellsize^2);
else
da1=max(Q.Z(ix1));
da2=max(Q.Z(ix2));
end
e1Ap=acos(((da1+da2)/da1)^-ref_theta);
e2Ap=acos(((da1+da2)/da2)^-ref_theta);
% Convert to degrees
e1Ap=rad2deg(e1Ap);
e2Ap=rad2deg(e2Ap);
% Howard 1971 slope method
ix3=S.IXgrid(ds{ii,1});
sl1=mean(G.Z(ix1));
sl2=mean(G.Z(ix2));
sl3=mean(G.Z(ix3));
r1=sl3/sl1;
r2=sl3/sl2;
if r1>1
e1Sp=0;
else
e1Sp=acos(r1);
end
if r2>1
e2Sp=0;
else
e2Sp=acos(r2);
end
% Convert to degrees
e1Sp=rad2deg(e1Sp);
e2Sp=rad2deg(e2Sp);
% Howard 1971 Minimum Power
if ~isempty(Q)
C1=da1*sl1;
C2=da2*sl2;
C3=(da1+da2)*sl3;
inner1=((C3^2)+(C1^2)-(C2^2))/(2*C1*C3);
inner2=((C3^2)+(C2^2)-(C1^2))/(2*C2*C3);
if inner1>1
e1Mp=0;
else
e1Mp=acos(inner1);
end
if inner>1
e2Mp=0;
else
e2Mp=acos(inner2);
end
% Convert to degrees
e1Mp=rad2deg(e1Mp);
e2Mp=rad2deg(e2Mp);
% Store out
pred_angle(ii,:)=[e1Ap e2Ap e1Sp e2Sp e1Mp e2Mp];
else
% Store out
pred_angle(ii,:)=[e1Ap e2Ap e1Sp e2Sp];
end
end
% Rotate all points to horizontal
[xu1r,yu1r]=rotcoord(xu1,yu1,-theta1ap,0,0);
[xu2r,yu2r]=rotcoord(xu2,yu2,-theta2ap,0,0);
[xdr,ydr]=rotcoord(xd,yd,-theta3ap,0,0);
% Simple approximation of linear segment on
% rotated positions
a1=xu1r\yu1r;
a2=xu2r\yu2r;
a3=xdr\ydr;
% Find projected y positions
yp1r=xu1r*a1;
yp2r=xu2r*a2;
yp3r=xdr*a3;
% Rotate back
[xp1,yp1]=rotcoord(xu1r,yp1r,theta1ap,0,0);
[xp2,yp2]=rotcoord(xu2r,yp2r,theta2ap,0,0);
[xp3,yp3]=rotcoord(xdr,yp3r,theta3ap,0,0);
% Calculate r-squared
r21=rsquared(yu1,yp1);
r22=rsquared(yu2,yp2);
r23=rsquared(yd,yp3);
% Convert to polar
[theta1,rho1]=cart2pol(xp1,yp1);
[theta2,rho2]=cart2pol(xp2,yp2);
[theta3,rho3]=cart2pol(xp3,yp3);
% Find maximum radii and thetas for those radii
[mrho1,ix1]=max(rho1);
[mrho2,ix2]=max(rho2);
[mrho3,ix3]=max(rho3);
theta1=theta1(ix1);
theta2=theta2(ix2);
theta3=theta3(ix3);
% Calculate interlink angles
[e1,e2,ba,split,d1,d2]=interangle(theta1,theta2,theta3);
e1=rad2deg(e1);
e2=rad2deg(e2);
ba=rad2deg(ba);
% Determine order of streams and define handedness
% sensu James & Krumbein, 1969
so1=max(so(us{ii,1}));
so2=max(so(us{ii,2}));
% so1 should be greater than or equal to so2 as output
% from nnodesupdown function
if so1>so2
if d1==1 & d2==2
hnd=1;
elseif d1==2 & d2==1
hnd=2;
elseif d1==1 & d2==1 & e1>e2
hnd=1;
elseif d1==1 & d2==1 & e1<e2
hnd=2;
elseif d1==2 & d2==2 & e1>e2
hnd=2;
elseif d1==2 & d2==2 & e1<e2
hnd=1;
else
hnd=3;
end
else
hnd=3;
end
% Convert direction angles to cardinal
card1=mod(-90-rad2deg(theta1),360);
card2=mod(-90-rad2deg(theta2),360);
card3=mod(-90-rad2deg(theta3),360);
link_angle(ii,:)=[S.x(con{ii}) S.y(con{ii}) ba e1 e2 split d1 d2 hnd so1 so2];
strm_dir(ii,:)=[card1 card2 card3];
fit_streams(ii,:)=[e1_dst numel(xu1) r21 e2_dst numel(xu2) r22 ds_dst numel(xd) r23];
else
link_angle(ii,:)=[S.x(con{ii}) S.y(con{ii}) NaN NaN NaN NaN NaN NaN NaN NaN NaN];
strm_dir(ii,:)=[NaN NaN NaN];
fit_streams(ii,:)=[NaN NaN NaN NaN NaN NaN NaN NaN NaN];
switch method
case {'area','slope','minimum_power'}
pred_angle(ii,:)=[NaN NaN];
case 'all'
if isempty(Q)
pred_angle(ii,:)=[NaN NaN NaN NaN];
else
pred_angle(ii,:)=[NaN NaN NaN NaN NaN NaN];
end
end
end
end
if make_shape
makejunctionshape(link_angle,pred_angle,strm_dir,fit_streams,method,1,out_dir,out_name);
end
junctions=makejunctiontable(link_angle,pred_angle,strm_dir,fit_streams,method);
else
num_n=numel(nlist);
% Preallocate cell array
junctions=cell(2,num_n);
for jj=1:num_n
nOI=nlist(jj);
% Preallocate arrays
la=zeros(num_con,11);
sd=zeros(num_con,3);
fs=zeros(num_con,9);
switch method
case {'area','slope','minimum_power'}
pa=zeros(num_con,2);
case 'all'
if isempty(Q)
pa=zeros(num_con,4);
else
pa=zeros(num_con,6);
end
end
for ii=1:num_con
if vld(ii)
% Extract ix lists and clip
conix=con{ii};
us1ix=us{ii,1};
us2ix=us{ii,2};
dsix=ds{ii,1};
if numel(us1ix)>nOI
us1ix=us1ix(1:nOI);
end
if numel(us2ix)>nOI
us2ix=us2ix(1:nOI);
end
if numel(dsix)>nOI
dsix=dsix(1:nOI);
end
% Extract x y node lists
% Confluence point
xc=xl(conix); yc=yl(conix);
% Upstream links
xu1=xl(us1ix); yu1=yl(us1ix);
xu2=xl(us2ix); yu2=yl(us2ix);
% Downstream link
xd=xl(dsix); yd=yl(dsix);
% Translate so all links originate from confluence
xu1=xu1-xc; xu2=xu2-xc;
yu1=yu1-yc; yu2=yu2-yc;
xd=xd-xc; yd=yd-yc;
%% Depending on orientation of links, doing a linear fit
% on the link may produce spurious results. To improve
% the result, each link is rotated towards horizontal,
% using the mean orientation of the link to find the angle
% of rotation
% Find mean angle of stream points from stream orientation
theta1ap=median(atan2(yu1,xu1));
theta2ap=median(atan2(yu2,xu2));
theta3ap=median(atan2(yd,xd));
% Extract and calculate stream distances
e1_dst=max(dst(us1ix))-dst(conix);
e2_dst=max(dst(us2ix))-dst(conix);
ds_dst=dst(conix)-min(dst(dsix));
% Calculate predicted angles
switch method
case 'area'
% Howard 1971 area method
ix1=S.IXgrid(us1ix);
ix2=S.IXgrid(us2ix);
if isempty(Q)
da1=max(A.Z(ix1).*A.cellsize^2);
da2=max(A.Z(ix2).*A.cellsize^2);
else
da1=max(Q.Z(ix1));
da2=max(Q.Z(ix2));
end
e1p=acos(((da1+da2)/da1)^-ref_theta);
e2p=acos(((da1+da2)/da2)^-ref_theta);
% Convert to degrees
e1p=rad2deg(e1p);
e2p=rad2deg(e2p);
% Store out
pa(ii,:)=[e1p e2p];
case 'slope'
% Howard 1971 slope method
ix1=S.IXgrid(us1ix);
ix2=S.IXgrid(us2ix);
ix3=S.IXgrid(dsix);
sl1=mean(G.Z(ix1));
sl2=mean(G.Z(ix2));
sl3=mean(G.Z(ix3));
r1=sl3/sl1;
r2=sl3/sl2;
if r1>1
e1p=0;
else
e1p=acos(r1);
end
if r2>1
e2p=0;
else
e2p=acos(r2);
end
% Convert to degrees
e1p=rad2deg(e1p);
e2p=rad2deg(e2p);
% Store out
pa(ii,:)=[e1p e2p];
case 'minimum_power'
% Howard 1971 minimum power method
ix1=S.IXgrid(us1ix);
ix2=S.IXgrid(us2ix);
ix3=S.IXgrid(dsix);
sl1=mean(G.Z(ix1));
sl2=mean(G.Z(ix2));
sl3=mean(G.Z(ix3));
q1=max(Q.Z(ix1));
q2=max(Q.Z(ix2));
q3=q1+q2;
C1=q1*sl1;
C2=q2*sl2;
C3=q3*sl3;
inner1=((C3^2)+(C1^2)-(C2^2))/(2*C1*C3);
inner2=((C3^2)+(C2^2)-(C1^2))/(2*C2*C3);
if inner1>1
e1p=0;
else
e1p=acos(inner1);
end
if inner>1
e2p=0;
else
e2p=acos(inner2);
end
% Convert to degrees
e1p=rad2deg(e1p);
e2p=rad2deg(e2p);
% Store out
pa(ii,:)=[e1p e2p];
case 'all'
% Howard 1971 area method
ix1=S.IXgrid(us1ix);
ix2=S.IXgrid(us2ix);
if isempty(Q)
da1=max(A.Z(ix1).*A.cellsize^2);
da2=max(A.Z(ix2).*A.cellsize^2);
else
da1=max(Q.Z(ix1));
da2=max(Q.Z(ix2));
end
e1Ap=acos(((da1+da2)/da1)^-ref_theta);
e2Ap=acos(((da1+da2)/da2)^-ref_theta);
% Convert to degrees
e1Ap=rad2deg(e1Ap);
e2Ap=rad2deg(e2Ap);
% Howard 1971 slope method
ix3=S.IXgrid(dsix);
sl1=mean(G.Z(ix1));
sl2=mean(G.Z(ix2));
sl3=mean(G.Z(ix3));
r1=sl3/sl1;
r2=sl3/sl2;
if r1>1
e1Sp=0;
else
e1Sp=acos(r1);
end
if r2>1
e2Sp=0;
else
e2Sp=acos(r2);
end
% Convert to degrees
e1Sp=rad2deg(e1Sp);
e2Sp=rad2deg(e2Sp);
% Howard 1971 Minimum Power
if ~isempty(Q)
C1=da1*sl1;
C2=da2*sl2;
C3=(da1+da2)*sl3;
inner1=((C3^2)+(C1^2)-(C2^2))/(2*C1*C3);
inner2=((C3^2)+(C2^2)-(C1^2))/(2*C2*C3);
if inner1>1
e1Mp=0;
else
e1Mp=acos(inner1);
end
if inner>1
e2Mp=0;
else
e2Mp=acos(inner2);
end
% Convert to degrees
e1Mp=rad2deg(e1Mp);
e2Mp=rad2deg(e2Mp);
% Store out
pa(ii,:)=[e1Ap e2Ap e1Sp e2Sp e1Mp e2Mp];
else
% Store out
pa(ii,:)=[e1Ap e2Ap e1Sp e2Sp];
end
end
% Rotate all points to horizontal
[xu1r,yu1r]=rotcoord(xu1,yu1,-theta1ap,0,0);
[xu2r,yu2r]=rotcoord(xu2,yu2,-theta2ap,0,0);
[xdr,ydr]=rotcoord(xd,yd,-theta3ap,0,0);
% Simple approximation of linear segment on
% rotated positions
a1=xu1r\yu1r;
a2=xu2r\yu2r;
a3=xdr\ydr;
% Find projected y positions
yp1r=xu1r*a1;
yp2r=xu2r*a2;
yp3r=xdr*a3;
% Rotate back
[xp1,yp1]=rotcoord(xu1r,yp1r,theta1ap,0,0);
[xp2,yp2]=rotcoord(xu2r,yp2r,theta2ap,0,0);
[xp3,yp3]=rotcoord(xdr,yp3r,theta3ap,0,0);
% Calculate r-squared
r21=rsquared(yu1,yp1);
r22=rsquared(yu2,yp2);
r23=rsquared(yd,yp3);
% Convert to polar
[theta1,rho1]=cart2pol(xp1,yp1);
[theta2,rho2]=cart2pol(xp2,yp2);
[theta3,rho3]=cart2pol(xp3,yp3);
% Find maximum radii and thetas for those radii
[mrho1,ix1]=max(rho1);
[mrho2,ix2]=max(rho2);
[mrho3,ix3]=max(rho3);
theta1=theta1(ix1);
theta2=theta2(ix2);
theta3=theta3(ix3);
% Calculate interlink angles
[e1,e2,ba,split,d1,d2]=interangle(theta1,theta2,theta3);
e1=rad2deg(e1);
e2=rad2deg(e2);
ba=rad2deg(ba);
% Determine order of streams and define handedness
% sensu James & Krumbein, 1969
so1=max(so(us1ix));
so2=max(so(us2ix));
% so1 should be greater than or equal to so2 as output
% from nnodesupdown function
if so1>so2
if d1==1 & d2==2
hnd=1;
elseif d1==2 & d2==1
hnd=2;
elseif d1==1 & d2==1 & e1>e2
hnd=1;
elseif d1==1 & d2==1 & e1<e2
hnd=2;
elseif d1==2 & d2==2 & e1>e2
hnd=2;
elseif d1==2 & d2==2 & e1<e2
hnd=1;
else
hnd=3;
end
else
hnd=3;
end
% Convert direction angles to cardinal
card1=mod(-90-rad2deg(theta1),360);
card2=mod(-90-rad2deg(theta2),360);
card3=mod(-90-rad2deg(theta3),360);
la(ii,:)=[S.x(conix) S.y(conix) ba e1 e2 split d1 d2 hnd so1 so2];
sd(ii,:)=[card1 card2 card3];
fs(ii,:)=[e1_dst numel(xu1) r21 e2_dst numel(xu2) r22 ds_dst numel(xd) r23];
else
la(ii,:)=[S.x(con{ii}) S.y(con{ii}) NaN NaN NaN NaN NaN NaN NaN NaN NaN];
sd(ii,:)=[NaN NaN NaN];
fs(ii,:)=[NaN NaN NaN NaN NaN NaN NaN NaN NaN];
switch method
case {'slope','area','minimum_power'}
pa(ii,:)=[NaN NaN];
case 'both'
if isempty(Q)
pa(ii,:)=[NaN NaN NaN NaN];
else
pa(ii,:)=[NaN NaN NaN NaN NaN NaN];
end
end
end
end
if make_shape
makejunctionshape(la,pa,sd,fs,method,jj,out_dir,out_name);
end
junctions{1,jj}=fit_distance(jj);
junctions{2,jj}=makejunctiontable(la,pa,sd,fs,method);
end % end nodes for
% Calculate sensitivity to fit distance
jamat=zeros(num_con,num_n);
e1mat=zeros(size(jamat));
e2mat=zeros(size(jamat));
for jj=1:num_n
J=junctions{2,jj};
jamat(:,jj)=J.junction_angle;
e1mat(:,jj)=J.e1_obs_angle;
e2mat(:,jj)=J.e2_obs_angle;
end
% Calculate mean and standard deviation
mean_junction_angle=mean(jamat,2);
mean_e1_angle=mean(e1mat,2);
mean_e2_angle=mean(e2mat,2);
std_junction_angle=std(jamat,1,2);
std_e1_angle=std(e1mat,1,2);
std_e2_angle=std(e2mat,1,2);
if strcmp(method,'slope') | strcmp(method,'all')