%Edwin Kite %******************* Import Data ************************** %1'XCoord' 2'YCoord' 3'FID' 4'ID' 5'PQUALITY' 6'PSTYLE' 7'TRAMBIG' %8'TRQUALITY' 9'SUREMEANDR' 10'ORIG_FID' 11'POINT_X' 12'POINT_Y' 13'038882HIZ' %14'548104HIZ' 15'B20G02CTXZ' a = dlmread('Transect_1_Channel_Centerlines_for_Third_Pass_Meander_Wavelength_Measurements_csv',',',1,0); %1 'XCoord' 2'YCoord' 3'FID' 4'ID' 5'ORIG_FID' 6'POINT_X' 7'POINT_Y' %8 '038882HIZ' 9'548104HIZ' 10'B20G02CTXZ a2 = dlmread('Transect_1_Pick_2_Third_Pass_Centerlines_for_Meander_Wavelength_Measurements_Points_csv',',',1,0); a3 = dlmread('Transect_1_Pick_3_Third_Pass_Centerlines_for_Meander_Wavelength_Measurements_csv',',',1,0); cid = a(:,10)+1; %ArcGIS is zero-indexed for feature IDs vx = a(:,1) - min(a(:,1)); %vertex x vy = a(:,2) - min(a(:,2)); %vertex y hiz = min(a(:,13),a(:,14)); ctxz = a(:,15); vz = min(ctxz,hiz); vz(vz==0) = NaN; dbip = 10; %distance between interpolated points (m) smoothn = 10; %number of interpolated points to smooth curvature for (to prevent artifactual inflection picks) minsinu = 1.1; %minimum half-meander sinuosity to be considered as a half-meander %********* Pepper Points Uniformly Along Channel Centerlines ************** ctl(1:max(cid)) = 0; %channel total length %Interpolate uniformly in s, dir coordinates for l = 1:(size(a,1)-1) if cid(l) == cid(l+1)%Is it from the same channel? lseg(l) = sqrt((vx(l+1) - vx(l)).^2 + (vy(l+1) - vy(l)).^2); %length of segment ctl(cid(l)) = ctl(cid(l)) + lseg(l); %channel total length lac(l+1) = ctl(cid(l)); %length along channel for this vertex else %last point on channel trace %lac(l) = ctl(cid(l)); end end for m = 1:max(cid) %position vector of interpolated points p(m).x = interp1(lac(cid==m),vx(cid==m),[0:dbip:ctl(m)]); p(m).y = interp1(lac(cid==m),vy(cid==m),[0:dbip:ctl(m)]); p(m).z = interp1(lac(cid==m),vz(cid==m),[0:dbip:ctl(m)]); p(m).dir = [ atan2d(p(m).y(2:end) - p(m).y(1:end-1),p(m).x(2:end) - p(m).x(1:end-1)) NaN]; %returns direction in degrees anticlockwise from E. Not defined for last point. p(m).th = [NaN p(m).dir(2:end) - p(m).dir(1:end-1) ]; %change in direction (degrees). Not defined for first or last points. p(m).th(p(m).th<-180) = 360 + p(m).th(p(m).th<-180); p(m).th(p(m).th> 180) = 360 - p(m).th(p(m).th> 180); p(m).smooth = smooth(p(m).th,smoothn)'; p(m).th1plusth2 = [NaN p(m).smooth(1:end-1) + p(m).smooth(2:end) ]; p(m).th2plusth3 = [ p(m).smooth(1:end-1) + p(m).smooth(2:end) NaN]; %Following Hemberger, 1986 Master's thesis, U. Virginia: p(m).isinfl = [1 - (sign(p(m).th1plusth2) == sign(p(m).th2plusth3))]; end scatter(vx,vy,100,'m*') %rem(a(:,10)*1e9,11) hold on for l = 1:size(p,2) scatter(p(l).x,p(l).y,100,p(l).isinfl,'Filled') end scatter([1:length(p(1).th)],p(1).th,100,p(1).isinfl,'Filled') hold on; title('debugging plot'); grid on plot(p(1).th) %******************* Define Half-Meanders ************************** for m = 1:max(cid) p(m).okhmep = p(m).isinfl; %O.K. half-meander end points p(m).okhmep(1) = 0; %no inflections at end of channel picks p(m).okhmep(end) = 0; %no inflections at end of channel picks hmeplist = find(p(m).okhmep==1); if length(hmeplist)>1 for epi = 1:(length(hmeplist)-1) %end-point index ep = hmeplist(epi); h(m).p(epi,1,1:2) = ... [p(m).x(hmeplist(epi)) p(m).x(hmeplist(epi+1))];%x,start|end h(m).p(epi,2,1:2) = ... [p(m).y(hmeplist(epi)) p(m).y(hmeplist(epi+1))];%y,start|end h(m).ls(epi) = ... %half-meander straight-line distance sqrt( ... (h(m).p(epi,2,2) - h(m).p(epi,2,1))^2 ... +(h(m).p(epi,1,2) - h(m).p(epi,1,1))^2); h(m).lc(epi) = (hmeplist(epi+1) - hmeplist(epi))*dbip;%along-curve distance h(m).th(epi) = nanmean(p(m).th(hmeplist(epi):hmeplist(epi+1))); h(m).sinu(epi) = h(m).lc(epi) ./ h(m).ls(epi); end %Half-meanders below a sinuosity threshold 'minsinu' will be rejected. %Remove end-points %that are the end-points for TWO half meanders below the sinuosity %threshhold from the end-point index %end-point index. %Then recalculate half-meander sinuosity between the remaining %end-points. %This done removing one end-point at a time in order of minimum %"average" sinuosity oksinu = h(m).sinu>minsinu; h(m).toberemoved = find(1+find((oksinu(1:end-1) + oksinu(2:end))<1)==1; [null min_mean_sinu] = min(h(m).sinu(1:end-1) + h(m).sinu(2:end)) if sum(min_mean_sinu==h(m).toberemoved) == 1 %minimum mean sinuosity is for an endpoint bracketed by two half-meanders with less than critical sinuosity? hmeplist(min_mean_sinu) = []; %half-meander end-point list end h(m).p = []; h(m).ls = []; h(m).lc = []; h(m).sinu = []; for epi = 1:(length(hmeplist)-1) %end-point index ep = hmeplist(epi); h(m).p(epi,1,1:2) = ... [p(m).x(hmeplist(epi)) p(m).x(hmeplist(epi+1))];%x,start|end h(m).p(epi,2,1:2) = ... [p(m).y(hmeplist(epi)) p(m).y(hmeplist(epi+1))];%y,start|end h(m).ls(epi) = ... %half-meander straight-line distance sqrt( ... (h(m).p(epi,2,2) - h(m).p(epi,2,1))^2 ... +(h(m).p(epi,1,2) - h(m).p(epi,1,1))^2); h(m).lc(epi) = (hmeplist(epi+1) - hmeplist(epi))*dbip;%along-curve distance h(m).sinu(epi) = h(m).lc(epi) ./ h(m).ls(epi); h(m).p(epi,3,1:2) = ... [p(m).z(hmeplist(epi)) p(m).z(hmeplist(epi+1))];%z,start|end end end end %****** Compute Summary Parameters & Plot Results for Debugging ********* figure mi = 0; mreji = 0; for m = 1:max(cid) %figure scatter(p(m).x,p(m).y,20,p(m).z,'Filled') hold on scatter(p(m).x(find(p(m).isinfl)),p(m).y(find(p(m).isinfl)),200,'k*') % scatter(p(m).x,p(m).y,100,p(m).isinfl,'Filled') hold on for n = 1:size(h(m).sinu,2) if h(m).sinu(n) > minsinu line([h(m).p(n,1,1) h(m).p(n,1,2)],[h(m).p(n,2,1) h(m).p(n,2,2)],'Color','g') text(h(m).p(n,1,1),h(m).p(n,2,1)+50,num2str(round(h(m).ls(n)))) mi = mi + 1; mls(mi) = h(m).ls(n); mx(mi) = mean(h(m).p(n,1,:)); my(mi) = mean(h(m).p(n,2,:)); mz(mi) = mean(h(m).p(n,3,:)); mm(mi) = m; msinu(mi) = h(m).sinu(n); else line([h(m).p(n,1,1) h(m).p(n,1,2)],[h(m).p(n,2,1) h(m).p(n,2,2)],'Color','r') mreji = mreji + 1; mrejsinu(mreji) = h(m).sinu(n); end title(num2str(m)) hold on axis equal end end figure scatter(mz,mls,30,mx,'Filled');grid on title('Half-meander straight-line distance versus elevation') figure hist([mrejsinu msinu],[1:0.01:4]) title('Histogram of all sinuosities')