 % This script uses stream vertices and samples riparian buffer cover along
 % the stream using transects running perpendicular to the stream.
 % requirements -- a streams file with all vertices IN ORDER spatially with
 % XY coordinates. An additional column specifies the width of the stream itself. That way the algorithm
 % doesn't start int he middle of a wide stream. For this study the width
 % was gathered from GIS methods given various water / landcover layers AND
 % manual interpretation / validation using high res imagery. 
 %The algorithm moves from vertext to vertex to measure
 % length along each stream segment. Ideally each stream has a unique id as
 % well. This will allow you to run particular subsets of the watershed.
 % you will also need to convert your discrete return lidar data into
 % matlab files for processing. Create two sets of files - one for the
 % vegetation and one for the ground points. The algorithm looks at the xy
 % extent of the stream reach being processed and cuts out the needed lidar
 % from the matlab matrices. This optimization allows you to process
 % streams more efficiently without looking at data beyond the reach study
 % extent (which is generally a 200 m buffer on either side of the stream
 % for this code). 
  
clc
clear
disp('BEGIN PROCESSING STREAM TRANSECTS')
tic

% ************************* Define All Variables *************************
conversion = 3.28083989501312; %conversion from meters to feet.
seg_length = 100; % distance between sample points along stream (m)

% distance between each sample point along perpendicular transect (m)
perp_seg = 1; % used to calculate buffer width 

transect_length = 200; % total distance to sample along transect (m)

% number of samples to take along perpendicular line (for buffer width calculation)
num_of_samples = floor((transect_length/perp_seg));

% sample square size (m) used to calculate the final lidar sample values % starting from the edge 
% of the buffer back to the stream
sampleSquare = 5; %this is the shorter edge - the thickness or height of the sample
sampleLength=10; %this is the length of the sample area (parallel to the stream)

% the distance between samples along the transect (added to allow for 10m sample box & 5m sampling rate
sampleDistance = 5;

% sampleBoxSize = window size used to INITIALLY calculate buffer width. 
% if perp_seg = 1, every 1 meter the script will sample lidar within a sampleBoxSize sized window
sampleBoxSize=5;
 
 % overstory height & fractional cover threshold [height,FC value]
 % These are the thresholds below which lidar ht and FC metrics are not considered "forest"
 thresholds=[6.56,.30];
 
 % 1= didn't run width; 0 = run buffer width
 dontRunWidth=0; 
 
 InteriorCellSize=10; % this is the cell size of the interior sample for the edge interior analysis

% ************************* READ IN STREAM DATA *************************

% Data should be formated in 5 columns: Name, VertexNum, X, Y, ChannelWidth)
% Name: 1 = Spring Creek, 2 = Tributary Spring Creek
% Stream data X,Y represents the XY location of each vertex along each stream centerline

fprintf('Import & Split Streams. \n\r')
%originalStream = dlmread('data\streamPoints_allClean4.csv',',',1,0); 
originalStream = dlmread('data\streamPoints_allClean4_shortList.csv',',',1,0); %a shorter list

% The code below will split stream data into smaller segments. If the stream segments are too long, 
% then there is too much to process efficiently. "Too much is dependent on computing resources.

% Added .1 m to this (split dist = 1000.1) to compensate for decimal effects resulting in
% 1 less transect on some reaches that are exactly 1000m long

%split into xm length segments - remove stream segments that are too short to sample
SplitStreams = funSplitStreams(1000.10, originalStream,seg_length); 

%split stream field 6 is the width of the stream in meters. 
save runOutput\SplitStreams SplitStreams
clear originalStream

% calculate length of new streams data along with number of vertices and total transects for that section of stream.
% AllStreamLength columns: StreamName, StreamLength(m), NumberOfVertices, Num of transects
AllStreamLength = funCalcStreamLength(seg_length, SplitStreams);
save runOutput\AllStreamLength AllStreamLength

%calculate total number of streams
totalStreams = length(AllStreamLength(:,1));

fprintf('Done reading in streams dataset!! \n\r')

%preallocate stream length matrix / calculate transect numbers per stream 
allTransects=sum(AllStreamLength(:,5));

%preallocate forest sample / calculate max number of transects
maxTransects=max(AllStreamLength(:,5));

%preallocate stream length / 
% 5 columns transectWidth=(:,:,j,Name)  %%DELETE THIS LINE??
StreamLength=zeros(totalStreams,2);

%for Name = 1:1
iteration=0;

% Use this to adjust which streams you wish to run in a particular iteration. this is useful 
% when you don't want to process an entire watershed. The code can be adjusted to run everything by 
% Setting startStream to 1 and commenting out totalStreams.
startStream=1; %this is  useful if you want to run a subset of streams in your dataset
%totalStreams=367;  

% Calculate total transects and streams for this run
y=0;
for u=startStream:totalStreams
    y=y+1;
    runTransects=(find(AllStreamLength(:,2)==u));
    transectNumber(y,:)=AllStreamLength(runTransects,5);
end

currentTotalTransects=sum(transectNumber);

%%


% ************************* BEGIN LOOP THROUGH STREAMS *************************

% this section processes each stream segment, breaks the stream into 100m
% segments and identifies the subsequent xy location of each transect.

% 7 columns / preallocate transect width matrix
TransectWidth= NaN((currentTotalTransects*2),7);


for Name = startStream:totalStreams
    
    %Process through (totalStreams) number of streams as defined / calculated above. 
    %Because the loop could start at stream 4, iteration keeps track of which loop iteration is being processed
     
    iteration=iteration+1;

    % **************** BEGIN calculate stream vertices and along stream sample pts ****************    
    fprintf('BEGIN Processing Next Stream...Stream-\n\r')
    disp(Name)

    % find vertices for each stream segment that makes up each stream  

    A = find(SplitStreams(:,2) == Name); 
    % populate original stream name for each segment
    originalStreamName=SplitStreams(A(1),1);

    % number of vertices in stream segment
    L = length(A);

    % Order of points
    Order = SplitStreams(A(1,1):A(end,1),3);

    % Find stream vertex coordinates
    Stream_temp = (SplitStreams(Order,4:6));% SplitStreams(Order,5)];

    % sort Stream Coord (not sure why this step is in here??)
    %Stream_coord = Stream_temp(Order,:);
    Stream_coord = Stream_temp;
       
    % preallocate stream segment for L (total stream vertices) and x and Y)
    p1_stream = zeros((L-1),3);
    p2_stream = zeros((L-1),3);
    dist_ft = zeros((L-1),1);
    dist_sum = zeros(L,1);
    
    % Solve for the distance between stream vertices
    fprintf('Calculating stream vertices for sampling...\n\r')

    for j = 1:L-1
        p1_stream(j,:) = (Stream_coord(j,1:3));   %first stream vertex coordinate
        p2_stream(j,:) = (Stream_coord(j+1,1:3)); %second stream vertex coordinate

        % distance between 2 points
        dist_ft(j) = sqrt((p2_stream(j,1)-p1_stream(j,1))^2+(p2_stream(j,2)-p1_stream(j,2))^2);

        % convert length from ft to meters
        dist = 0.3048*dist_ft;
    end

    % Calculate cummulative stream length starting at first point given in file
    for k = 2:L-1
        dist_sum(1) = 0;
        dist_sum(2) = dist(1); % pull the distance between the vertices (dist matrix)
        dist_sum(k+1) = dist_sum(k)+dist(k);
    end

    %record stream length
    StreamLength(Name,:)= [Name, (max(dist_sum)) ];

    %(floor = round down); calculate number of stream transects per reach
    Num_of_trans = AllStreamLength(Name,5);
    seg=zeros((Num_of_trans),1); %preallocate for seg matrix
    
    % locate distance where the transect should begin
    for n = 1:Num_of_trans;
        seg(n) = find(dist_sum < seg_length*n, 1, 'last' ); 
    end

   
    % pull out points bounding each stream segment end coordinate

    % **************** END calculate stream vertices and sample pts along stream ****************    

    % check to see if number of transects is less than 1
    if Num_of_trans < 1
        fprintf('Stream sample distance is too long for the length of stream!.\n')
    else
     
    %if output values are off it could be because of this preallocation clearing out variables it needed.    
    p1_seg=NaN((length(seg)),3); %preallocate for seg matrices
    p2_seg=NaN((length(seg)),3);
    seg_dist=NaN((length(seg)),1);
    alpha=NaN((length(seg)),1);
    m=NaN((length(seg)),1);
    b=NaN(Num_of_trans,1);
    y1=NaN(Num_of_trans,1000);
    y2=NaN(Num_of_trans,1000);
    bufferSamplePoly=NaN((num_of_samples*2),11,Num_of_trans);
    pn=NaN(Num_of_trans,6);

    %calculate transects for stream
    %Num_of_trans=3 %use to test if you want to just run a few transects on each stream 
    
    for i = 1:Num_of_trans; 
        p1_seg(i,:) = p1_stream(seg(i),:);
        p2_seg(i,:) = p2_stream(seg(i),:);

        % calculate the distance along the line between p1 and p2 that = seg_length
        seg_dist(i) = seg_length*i - dist_sum(seg(i));

        % calculate slope between p1 and p2 for seg(i) 
        % if the line is vertical (the x values are the same) then the slope can't be calculated
        if p2_seg(i,1) == p1_seg(i,1)
            m(i,1)=NaN;
        else
            m(i,1) = (p2_seg(i,2)-p1_seg(i,2))/(p2_seg(i,1)-p1_seg(i,1));
            % Define whether the line is moving north, south, east or west
            % 1= the line is moving either east or north, 0 is west or south

            %for x coordinates (east, west)
            if p2_seg(i,1) > p1_seg(i,1)
                m(i,2) = 1;
            else
                m(i,2) = 0;
            end
            
            %for y coordinates (north, south)
            if p2_seg(i,2) > p1_seg(i,2)
                m(i,3) = 1;
            else
                m(i,3) = 0;
            end
        end

        % convert slope to angle (radians)
        alpha(i) = atan(m(i,1));

        % calculate new XY values which will be the starting point for the transect
        % to be drawn at the location Xm from the previous transect (defined above)
        [px, py]=funCalculateNewXY(p1_seg(i,:), p2_seg(i,:), seg_dist(i), alpha(i));

        %calculate average stream width at px,py
        avWidth=( (p1_seg(i,3)+p2_seg(i,3)) /2 );

        % define a range of x
        x = linspace(0,1e7,1000);

        % calculate the y intercept of the line
        b(i) = p1_seg(i,2)-m(i,1)*p1_seg(i,1);

        % solve for y (line between stream verticies)
        y1_temp = m(i,1)*x+b(i);

            % solve for the line orthogonal to the stream
            % Fixed to account for m=0 which was bombing the code
            if m(i,1) == 0
                y2_temp=py;
            else 
                y2_temp = (-1/m(i,1))*(x-px)+py;
            end
            
        % Assemble all variables into matricies
        y1(i,:) = y1_temp; 
        y2(i,:) = y2_temp;

        % coordinate of point = seg_dist
        pn(i,:) = [px py m(i,:) avWidth]; %adding av width to this here
        Pn(i,:,Name) = pn(i,:);
    end

    %preallocate for pnperp
    pn_perp_plot=zeros((num_of_samples*2),4,Num_of_trans);
    %pn_perp_plot=zeros((num_of_samples*2),2,Num_of_trans,totalStreams);


    %Solve for points along line perpendicular to stream segment at point Pn
    fprintf('Running transects...\n\r')

    %preallocate for the slope matrix
    m_perp=zeros((length(seg)),3);


    % **************** CALCULATE BUFFER SAMPLE POINTS ALONG TRANSECTS ****************

    %for j = 1:length(seg)
    for j = 1:Num_of_trans

        % calculate the slope of line perpendicular to stream segment
        m_perp(j,:) = [(-1./pn(j,3)) pn(j,4) pn(j,5) ];

        % convert slope to angle (radians)
        % alpha_perp(j) = atan(m_perp(j)); %put this in the function as it doesn't need to be calculated here
        Sampling(Name,:,:) = [Name;(Num_of_trans);(num_of_samples*2);]; %might be able to move this to the top

        % find point coordinates along line "perp_seg" distance from Pn
        % below is how far out from the stream along the transect it samples
        startingPoint=[pn(j,1) pn(j,2)];
        widthAtTransect=(pn(j,6));

        % 8/4/2011 Modification: consider stream width when calculating
        % points... so go out from the stream xm before starting the transect...

        %calculate average stream width between the two vertices

        [pn_perp_pos, pn_perp_neg] = funCalculatePointsAlongLine(startingPoint,num_of_samples,m_perp(j,:),perp_seg,widthAtTransect);

        % Assemble sample points north and south of stream into one matrix
        % 1:positive, 0: negative (below stream)
        %PN Perp_Plot Columns: Pos/Neg(which side of stream); ID, X, Y -- i think this is the
        %these xy locations are saved in the forestSample matrix below...
        pn_perp_plot(:,:,j) = [pn_perp_pos; pn_perp_neg];

        %for s = 1:pn_perp_length
        totalSamples2=(num_of_samples*2) ;
        bufferSamplePoly_temp2=zeros(totalSamples2,11);

    %**************** CALCULATE Sample Point Polygons for transect sample locations (j) ****************

        for s= 1:totalSamples2
            % calculate xm (samplesizebox) square vertices that will be used to sample forest along the
            samplePtBufWidth_temp=[pn_perp_plot(s,3,j) pn_perp_plot(s,4,j)];

            %calculate box vertex locations (where the lidar will be sampled)
           % bufferSamplePoly_temp = calculateBoxVertices(m_perp(j), samplePtBufWidth_temp, sampleBoxSize);
            %changed this to also use a rectangular sample area
            bufferSamplePoly_temp = funCalculateRectangleVertices(m_perp(j), samplePtBufWidth_temp, sampleBoxSize, sampleLength); 


           % add the stream name, transect number and direction (up or down along
           % the stream to and populate a matrix for the entire transect
            bufferSamplePoly_temp2(s,:)=[Name j pn_perp_plot(s,1,j) bufferSamplePoly_temp]; %might not need Name and j here...   
        end

        %this is the final matrix that contains all of the sample box vertex
        %points for the current transect (j) 
        bufferSamplePoly(:,:,j)=(bufferSamplePoly_temp2);


    end
        %save out bufferSamplePoly for future use
        %Columns: StreamName, transectNumber, Pos/Neg , then the XY coordinate locations for
        %the associated sample box
        saveBufSam=(['save intResults/bufferSamplePoly' num2str(startStream) '_' num2str(totalStreams) ' bufferSamplePoly']);
        eval(saveBufSam);

%%
    % **************** LOAD LIDAR DATA For STREAM SEGMENT BEING PROCESSED ****************
    %load lidar data needed for processing based upon the extent of the streams data and calculate value within each 1x1m cell 
    %LiDAR data should be processed and converted into XYZ format. For this code, save the vegetation points in one file
    % and the ground points in another file 
    % You can export data to XYZ with tools like lastools and many other tools. The data were
    %then imported into matlab and exported as .mat files to expedit load
    %time.
    
    %import lidar data to be used
    disp('LOADING LIDAR to calculate buffer width')
    tic
    disp('Now finding veg lidar extent and then populating veg matrix...')
    streamVegLidar=funFindStreamVegLidarExtent(Stream_temp);
    fprintf('Veg Lidar Data Processing complete, memory cleared...\n\r')

    toc
    tic
    disp('Now finding GROUND lidar extent and then populating GROUND matrix... ')
    streamGrdLidar=funFindStreamGrdLidarExtent(Stream_temp);
    fprintf('Ground Lidar Data Processing complete, memory cleared...\n\r')
    toc

    disp('Begin CALCULATE buffer width for each transect')
    fprintf('Find the buffer edge - sample lidar part I...\n\r')

    % calculate for each point along the transect whether it is forest or non
    % forest


    %preallocate ForestSample & sample temp
    %forestSample_temp is the matrix that contains the xy values of each sample
    %location and the associated value (forest vs non forest) calculated
    forestSample_temp=NaN((num_of_samples*2),7);
    forestSample=NaN((num_of_samples*2),7,Num_of_trans);
    for j=1:Num_of_trans 

        disp('transect number')
        disp(j)
        tic
        [forestSample_temp, veg, grd]=funSampleLidar_forestClass(bufferSamplePoly(:,:,j),streamVegLidar,streamGrdLidar,pn_perp_plot(:,:,j));
        forestSample(:,:,j)=(forestSample_temp); 
        toc 
    end

        % save out the forestSampleData
        %Columns = StreamName; TransectNumber; Pos or Neg (0 1) - which side of the stream;XCoord; YCoord; Ht (70th pct?),
        %FractionalCover
        saveForestSample=(['save results/forestSample_Strm' num2str(Name) ' forestSample']);
        eval(saveForestSample);

    % -----------------------------------------
    startingPoint2=[Pn(j,1,Name) Pn(j,2,Name)];

    % -------------------------CALCULATE BUFFER WIDTH

    % pn_perp_length = length(pn_perp); %not using this anywhere...
    %% 

    fprintf('Calculate buffer width (along transect)...\n\r')

    %------------ CALCULATE BUFFER WIDTH POINTS ALONG TRANSECTS -------------

    transectWidth_temp=NaN((Num_of_trans*2),14);
    l=1;

    %clear out lidar sample location
    lidarSampleLoc=[];

    for j = 1:Num_of_trans

       % disp('transect number ')
       % disp(j)

    % -----------------------------------------
    startingPoint2=[Pn(j,1,Name) Pn(j,2,Name)];
    currentTransWidth2=(Pn(j,6,Name));
    %}

    tempLength=length(forestSample(:,:,j));
    halfTempLength=(tempLength/2);

    %forestSample(:,:,j)
    BufferSample_pos = (forestSample(1:halfTempLength,:,j));
    BufferSample_neg = (forestSample((halfTempLength+1):end,:,j)); 


    %calculate  forest buffer width for each transect north and south of the stream
    %added transect begin xy at the end of this ... (2 more columns) Added
    %StreamWidth (currentTransWidth2); slope[m_perp(j)] 
    
    transPos_temp = [originalStreamName Name j 1 (funCalcTransectBufferWidth2(BufferSample_pos, perp_seg, startingPoint2,thresholds)) currentTransWidth2 m_perp(j)];
    transNeg_temp = [originalStreamName Name j 0 (funCalcTransectBufferWidth2(BufferSample_neg, perp_seg, startingPoint2,thresholds)) currentTransWidth2 m_perp(j)];


    %*********************look at this - i think the save should go at the end of this loop...
    transectWidth_temp(l,:)=(transPos_temp);
    transectWidth_temp(l+1,:)=(transNeg_temp);
    l=l+2;
    %this matrix contains the beginning and the end of the stream x and y value
    %columns: OriginalStreamName, Stream Name, TransectNum, pos/neg, Xbegin,
    %YBegin, Xend, Yend, finalWidth, density for the width span (FC along
    %transect), XY (starting point??) <will add m_perp(j) to this to have slope
    %as well>>
    saveTran=(['save results/transectWidth' num2str(Name) ' transectWidth_temp']);

    eval(saveTran);
    %m_perp(j)
    saveSlope=(['save results/slope' num2str(Name) ' m_perp']);
    eval(saveSlope);
    %%

    % this will load the slope values for the current stream segment if i
    % didn't already calculate stream width
    if dontRunWidth == 1
    loadSlope=(['load results/slope' num2str(Name) ' m_perp']);
    eval(loadSlope);
    end

    %calculate the total number of sample boxes that will fit within the buffer
    %(round down)

    %the edge of the buffer, x, y and the total buffer width (if it equals 0,
    %don't calculate lidar values)
    bufferEdgeLocPos=(transPos_temp(7:9));
    bufferEdgeLocNeg=(transNeg_temp(7:9));

    %make sure to include code of length = 0 or is too small to sample...
    %calculate the total number of samples for the transect 
    %am going to subtract a small amount to account for the edge effects at the
    %buffer edge (where the canopy may overlay the fields)
    edgeEffect=3;
    % calculate number of samples by subtracting edge effects from the buffer
    % width (transPos_temp(9)), then dividing by the square sample
    if sampleDistance==sampleSquare
        lidarSamplingPos=floor(((transPos_temp(9) - edgeEffect))/sampleDistance);
        lidarSamplingNeg=floor(((transNeg_temp(9) - edgeEffect))/sampleDistance);
    else
        difference=((sampleSquare)/2);
        lidarSamplingPos=floor(((transPos_temp(9) - edgeEffect-difference))/sampleDistance);
        lidarSamplingNeg=floor(((transNeg_temp(9) - edgeEffect-difference))/sampleDistance);  
    end

    %calculate transect sample center point
    transCenterDistPos=(transPos_temp(9) / 2);
    transCenterDistNeg=(transNeg_temp(9) / 2);

    %calculate samplingLocations from edge down to stream
    %num of sampling points is lidarSamplingPos or neg 

    %calculatePointsAlongLine (startPointPos, startPointNeg, sampleNumPos, sampleNumNeg,slope, samplingDistance)
    %disp('Calculating  box to run the final lidar sample along the transect for the width of the buffer (calculatePointsAlongLine_Reverse)')
    disp(j) % display current transect number
    %clear out lidarSample Loc
    lidarSampleLoc_temp=[];
    lidarSampleLocInterior_temp=[];
    %calculate all points along transect to sample
    %%
    %outputs are Transect Name Name; pos neg; buffer width; distance from
    %center and xy location of sample point
    %october 31 - changed sampleSquare to sample distance below -- that way i
    %get locations along the transect every "sample distance" distance however
    %i can make the sample area larger.
    %edgeEffect=0; %changing this oct 31 given the wider square size
    [lidarSampleLoc_temp] = funCalculatePointsAlongLine_Reverse(bufferEdgeLocPos, bufferEdgeLocNeg,lidarSamplingPos,lidarSamplingNeg,m_perp(j),sampleDistance,Name,j,edgeEffect,sampleSquare);

    %calculate interior / ext/center sample point locations
        %1 = interior, 2=bufferedge (landuse) 3= stream edge(stream) 
        %4 = Matrix Fringe, 5 stream fringe, 

    %if  buffer width is less than XXX then don't calculate points??
    [lidarSampleLocInterior_temp] = funCalculatePointsAlongLine_Reverse_TransectCenter(bufferEdgeLocPos, bufferEdgeLocNeg,1,1,m_perp(j),transCenterDistPos,transCenterDistNeg,Name,j,edgeEffect,sampleSquare);


    % calculate box vertices for interior transect location
    for h=1:length(lidarSampleLocInterior_temp(:,1)) % This should always = 10 changed to 10 nov 1 when i added "fringe" values
        %as 5 samples per transect - 2 transects (1 in each direction) to only run samples for sample locations that do not equal the
        %starting point...
        if h==1 || h==6
            cellSize=InteriorCellSize;
        else
            cellSize=sampleSquare;
        end
        %only draw a box if the buffer width is larger than 25m - otherwise
        %there is not enough space to sample.
        %if lidarSampleLocInterior_temp(h,4) > 20

        sample_point=(lidarSampleLocInterior_temp(h,6:7));     
        %sampleInteriorPolyVert_temp2 = cal1culateBoxVertices(m_perp(j), sample_point, cellSize);
        sampleInteriorPolyVert_temp2 = funCalculateRectangleVertices(m_perp(j), sample_point, cellSize, sampleLength); %changed to allow for a rectangular selection
        %end
        %add stream name, transect, north/south etc to matrix
        %Name; TransectNum (J); Pos/Neg ; BufferWidth; (xy location of center sample); DistFromStream (1/2 buffer width); XY locations of the 4 vertices of the sample polygon; 
        sampleInteriorPolyVert_temp(h,:)= [Name j lidarSampleLocInterior_temp(h,3:4) lidarSampleLocInterior_temp(h,6:7) lidarSampleLocInterior_temp(h,5) sampleInteriorPolyVert_temp2];
    end
 
    %adjusted the code above to buffer the sampling starting at the stream edge
    %(measured in arcmap using imagery). it is a rough estimate but will bring
    % also added an edge effect to consider tree overhang buffer width to within a 1-3m accuracy where the stream is wide.
    % Also need to account for samples where
    %width = 0 (make sure it's not a part of the lidar sampling.

    %for all points, calculate
    sizeSam=length(lidarSampleLoc_temp(:,1));
    samplePolyVert_temp=zeros(sizeSam,15);

    for h=1:length(lidarSampleLoc_temp(:,1))
        %need to only run samples for sample locations that do not equal the
        %starting point...
        sample_point=(lidarSampleLoc_temp(h,6:7));


        %sampleLength = the length of the sample - i will be sampling parallel to
        %the stream xm (10-20m) by 5m tall; sampleSquare is the height or depth
        %of each sample area (5m)

        samplePolyVert_temp2 = funCalculateRectangleVertices(m_perp(j), sample_point, sampleSquare, sampleLength); 
        %add stream name, transect, north/south, sample XY (center point (lidarsampleLoc 6:7), and distance measured (lidarSample(5)  etc to matrix
        samplePolyVert_temp(h,:)= [Name j lidarSampleLoc_temp(h,3:4) lidarSampleLoc_temp(h,6:7) lidarSampleLoc_temp(h,5) samplePolyVert_temp2];
    end

        %clear out lidar sample location (the coordinates of the sample point
        %running from the buffer edge back towards the stream)

    if iteration == 1 && j == 1  %for the first transect / stream

            sampleInteriorPolyVert=(sampleInteriorPolyVert_temp);
            samplePolyVert=(samplePolyVert_temp);
            %Columns: streamName; TransectNumber; Pos/Neg; BufferWidth; Distance from Edge; X, Y
            lidarSampleLoc=(lidarSampleLoc_temp);
            % lidarInteriorSampleLoc=(lidarInteriorSampleLoc_temp);

        else
            lidarSampleLoc=[lidarSampleLoc;lidarSampleLoc_temp];      
            samplePolyVert=[samplePolyVert;samplePolyVert_temp];
            sampleInteriorPolyVert=[sampleInteriorPolyVert;sampleInteriorPolyVert_temp];

    end
        %save out results of lidar sampling vertices
        %these are the sample boxes for the final transect that begin at
        %width(x) calculated previously and continue all the way to the stream
        %center point...


    %%
    end 


    end
    %%
    % extract polygon vertices for the current stream (name = x) both for all
    % transect samples and for the interior sample.
    % 
    %first calculate for all transect samples

        samplePolyVert_StrmIndex = (find(samplePolyVert(:,1)==Name)) ;
        sampleNumbers=length(samplePolyVert_StrmIndex);
        samplePolyVert_StrmCurrent=zeros(sampleNumbers,15);

        for polyIndex=1:sampleNumbers
            samplePolyVert_StrmCurrent(polyIndex,:) = (samplePolyVert(samplePolyVert_StrmIndex(polyIndex),:)); 
        end

        %then calculate values within polygon for the single interior polygon

        sampleInteriorPolyVert_StrmIndex = (find(sampleInteriorPolyVert(:,1)==Name)) ;
        sampleInteriorNumbers=length(sampleInteriorPolyVert_StrmIndex);

        sampleInteriorPolyVert_StrmCurrent=zeros(sampleInteriorNumbers,15);
        for polyInteriorIndex=1:sampleInteriorNumbers
            sampleInteriorPolyVert_StrmCurrent(polyInteriorIndex,:) = (sampleInteriorPolyVert(sampleInteriorPolyVert_StrmIndex(polyInteriorIndex),:)); 
        end

       %save polygon vertices that are used to define the sample area to calculate final lidar values
       % within the x by xm squares in case i want to use them later... 
       savePolyVert=(['save results/samplePolyVert' num2str(Name) ' samplePolyVert_StrmCurrent']);
       eval(savePolyVert); 

       saveInteriorPolyVert=(['save results/sampleInteriorPolyVert' num2str(Name) ' sampleInteriorPolyVert_StrmCurrent']);
       eval(saveInteriorPolyVert); 

    %%
    % begin sample lidar data within interior and all polygons along transect...

    % check the lidar sample location matrix - i may no longer need it...

    % sample ALL polygons first

    % Make sure lidar sample results temp is cleared out
    lidarSampleResults_temp=[];
    % sample lidar data within all polygon sample boxes (for the final width of the stream, and for all transects in the current stream reach
    %the last number is the width threshold (the value below which values will
    %not be calculated

    [lidarSampleResults_temp vegLidar_temp grdLidar_temp]= funSampleLidar(samplePolyVert_StrmCurrent, streamVegLidar, streamGrdLidar, 8);
    disp('Done calculating sample lidar values for transect')

    %sample interior polygons second
    %lidarSampleLoc is just the XY value for the stream center <this needs to be fixed>
    % Make sure lidar sample results temp is cleared out
    lidarInteriorSampleResults_temp=[];
    % sample lidar data within all polygon sample boxes (for the final width of the stream, and for all transects in the current stream reach
    [lidarInteriorSampleResults_temp vegInteriorLidar_temp grdInteriorLidar_temp]= funSampleLidar(sampleInteriorPolyVert_StrmCurrent, streamVegLidar, streamGrdLidar,25);

    disp('Done calculating sample lidar values for interior / edge for transect')

    %write out sample results to a matrix
        if iteration == 1
           % VegSampleLidar=(vegLidar_temp); %not going to save this each time as it will really get large
           % GrdSampleLidar=(grdLidar_temp);
            LidarSampleResults= (lidarSampleResults_temp);
            LidarInteriorSampleResults= (lidarInteriorSampleResults_temp);
        else
            LidarSampleResults = [LidarSampleResults;lidarSampleResults_temp]; 
            LidarInteriorSampleResults = [LidarInteriorSampleResults;lidarInteriorSampleResults_temp]; 
           % VegSampleLidar=[VegSampleLidar;vegLidar_temp];
           % GrdSampleLidar=[GrdSampleLidar;grdLidar_temp];    
        end

        % save out lidar results 
        saveVegLidar=(['save results/lidar/vegLidar_Strm' num2str(Name) ' vegLidar_temp']);
        eval(saveVegLidar);
        saveGrdLidar=(['save results/lidar/grdLidar_Strm' num2str(Name) ' grdLidar_temp']);
        eval(saveGrdLidar);

        %save out interior edge lidar
        saveVegLidar=(['save results/lidar/vegInteriorLidar_Strm' num2str(Name) ' vegInteriorLidar_temp']);
        eval(saveVegLidar);
        saveGrdLidar=(['save results/lidar/grdInteriorLidar_Strm' num2str(Name) ' grdInteriorLidar_temp']);
        eval(saveGrdLidar);

        saveLidarResults=(['save results/lidarSampleResults_Strm' num2str(Name) ' lidarSampleResults_temp']);
        eval(saveLidarResults);
        saveLidarInteriorResults=(['save results/lidarInteriorSampleResults_Strm' num2str(Name) ' lidarInteriorSampleResults_temp']);
        eval(saveLidarInteriorResults);

        %going to keep this temporarily ...
         %clear streamGrdLidar
         %clear streamVegLidar



        %%
        %  
        if iteration ==1
          TransectWidth=(transectWidth_temp);
        else
          TransectWidth=[TransectWidth;transectWidth_temp];
        end
     %%
 end
 
%%
%write out data to be plotted and exported to text
%calculate variables to preallocate matrices
tran=1;
%it's ok to use Sampling here instead of All stream length because Sampling
%just calculates values for what stream segments were run (whereas
%allStreamLength calculates for ALL streams)
samplesPerTransect=((num_of_samples*2)); %total number of points per transect (for both positive and negative along the transect)
totalTransects=sum(Sampling(:,2)); % total number of transects in the entire dataset %might be able to use allTransects already calculated
totalPolygonPoints= ((totalTransects*samplesPerTransect)*4);
totalStreams = length(Sampling(:,1)); %total number of streams
totalSamples= totalTransects * (Sampling(:,3)); %total number of samples taken throuhout the watershed

g=1;
p=1;
l=1;

disp('done writing out data')
%%
fprintf('Done Processing streams:\n\r')
disp(startStream)
disp(totalStreams)

%%


% ***************** write out final Results transect values *****************


%---------------------------turn off plot begin

figure(3)

%Move this into a plot function for later
 plot(vegLidar_temp(:,1),vegLidar_temp(:,2),'ro','MarkerFaceColor','r','MarkerSize',3), hold on
% plot(grdLidar_temp(:,1),grdLidar_temp(:,2),'k*'), hold on

%set('MarkerEdgeColor','r','MarkerFaceColor','none') %use this to fill in
%symbols in a plot
%{
% *********************Old Sample Data to Plot
plot(SplitStreams(1:10,5),SplitStreams(1:10,6),'--bo','LineWidth',1.25), hold on
plot(SplitStreams(11:24,5),SplitStreams(11:24,6),'--bo','LineWidth',1.25), hold on
% *********************************************
%}
%
%
    plot(samplePolyVert(5,8),samplePolyVert(5,9),'r*','MarkerFaceColor','r'),hold on
    plot(samplePolyVert(5,10),samplePolyVert(5,11),'r*','MarkerFaceColor','r'),hold on
    plot(samplePolyVert(5,12),samplePolyVert(5,13),'r*','MarkerFaceColor','r'),hold on
    plot(samplePolyVert(5,14),samplePolyVert(5,15),'r*','MarkerFaceColor','r'),hold on
%}
%{
% Plot samples (not edges)
polyNoEdge=(samplePolyVert((samplePolyVert(:,7)~=999),:));
for i=1:length(polyNoEdge)
    plot(polyNoEdge(i,8),polyNoEdge(i,9),'ko','MarkerFaceColor','k'),hold on
    plot(polyNoEdge(i,10),polyNoEdge(i,11),'ko','MarkerFaceColor','k'),hold on
    plot(polyNoEdge(i,12),polyNoEdge(i,13),'ko','MarkerFaceColor','k'),hold on
    plot(polyNoEdge(i,14),polyNoEdge(i,15),'ko','MarkerFaceColor','k'),hold on
end

%plot samples (edges only)
polyEdge=(samplePolyVert((samplePolyVert(:,7)==999),:));
for i=1:length(polyEdge)
    plot(polyEdge(i,8),polyEdge(i,9),'bo','MarkerFaceColor','y'),hold on
    plot(polyEdge(i,10),polyEdge(i,11),'bo','MarkerFaceColor','y'),hold on
    plot(polyEdge(i,12),polyEdge(i,13),'bo','MarkerFaceColor','y'),hold on
    plot(polyEdge(i,14),polyEdge(i,15),'bo','MarkerFaceColor','y'),hold on
end
    %plot sample locations
    plot(LidarSampleResults(:,5), LidarSampleResults(:,6),'bs','MarkerSize',5), hold on
% plot sample locations
%}


%
plot(sampleInteriorPolyVert(:,8),sampleInteriorPolyVert(:,9),'bo','MarkerFaceColor','r'),hold on
plot(sampleInteriorPolyVert(:,10),sampleInteriorPolyVert(:,11),'bo','MarkerFaceColor','r'),hold on
plot(sampleInteriorPolyVert(:,12),sampleInteriorPolyVert(:,13),'bo','MarkerFaceColor','r'),hold on
plot(sampleInteriorPolyVert(:,14),sampleInteriorPolyVert(:,15),'bo','MarkerFaceColor','r'),hold on
%}
%plot the center of the transect
plot(sampleInteriorPolyVert(:,5),sampleInteriorPolyVert(:,6),'b*','MarkerFaceColor','r'),hold on


plot(TransectWidth(:,7), TransectWidth(:,8),'m*','MarkerSize',30), hold on
%instead of this plot using the LidarSampleResults matrix


allStreamsLen=(length(AllStreamLength(:,1)));


for f=1:allStreamsLen
 %for f=8:8
    vertices=AllStreamLength(f,4);
    A2=zeros(vertices);
    A2 = find(SplitStreams(:,2) == f); 
    % number of verticies in stream segment
    L2 = length(A2);
    plot(SplitStreams((A2(:,1)),4),SplitStreams((A2(:,1)),5),'--bo','LineWidth',1.25,'MarkerSize',5), hold on
end

%{
%plot sample locations going away from the stream (used initially to
%determine buffer width
for streamLoop=1:totalStreams  
 for chaya=1:Sampling(streamLoop,2)
 plot(pn_perp_plot(:,3,chaya,streamLoop), pn_perp_plot(:,4,chaya,streamLoop),'ks'), hold on
 line(pn_perp_plot(:,3,chaya,streamLoop),pn_perp_plot(:,4,chaya,streamLoop)), hold on
 end
end 
%}

%{
% plot DEM sample selected (closest point locations)
for streamLoop=1:totalStreams
 for chaya2=1:Sampling(streamLoop,2) %for each transect
 plot(CHMSample(:,1,1,chaya2,streamLoop), CHMSample(:,2,1,chaya2,streamLoop),'go'), hold on
 end
end 
%}

%{
for streamLoop=1:totalStreams
 for chaya3=1:Sampling(streamLoop,2) %for each transect
 plot(PolygonVertices(:,1,1,chaya3,streamLoop), PolygonVertices(:,2,1,chaya3,streamLoop),'b*'), hold on
 plot(PolygonVertices(:,3,1,chaya3,streamLoop), PolygonVertices(:,4,1,chaya3,streamLoop),'b*'), hold on
 plot(PolygonVertices(:,5,1,chaya3,streamLoop), PolygonVertices(:,6,1,chaya3,streamLoop),'b*'), hold on
 plot(PolygonVertices(:,7,1,chaya3,streamLoop), PolygonVertices(:,8,1,chaya3,streamLoop),'b*'), hold on
 end
end 
%}

%transectWidth(:,:,j,Name)
axis equal
%axis([ 1950000 1954000 249000 253000])
set(gca,'YTickLabel',num2str(get(gca,'YTick').'));
set(gca,'XTickLabel',num2str(get(gca,'XTick').'));
     

