% This code will plot out streams data for as many or as few streams as
% indicated by the start and end stream variables

%the output is as follows
%ForestSample: 
    %this matrix contains the xm sample point locations (xy) and
    %associate FC and av height value for the entire length of each transect.
    %it is a multi dimensional matrix with each transect being a difference
    %dimension (400, 7 , 9) = 400 samples along the transect (both directions,
    %in total) 7 columns, 9 transect in total)
    
%the finalTransects output includes the stats to be used for final analysis
%
clc
clear


startStream=1;
endStream=1; %this number will vary depending on how many streams are being processed
%{
startStream=1;
endStream=367;
%}
sampleXY=[];
InteriorEdge=[];
TransectSamples_UniqueID=[];
FinalTransects=[];

load runOutput/SplitStreams
load runOutput/AllStreamLength

% this file was derived using ArcGIS methods. Each transect was classified
% based upon the landuse type adjacent to it. We used local landuse data
% for this step as it was the more current and accurate layer available and
% we had data for the same time the lidar data were flown. THe matrix
% should be an X (number of transects) x 2 dimension matrix. where column
% one is the matrix ID and column 2 is the landuse class. Landuse classes
% for this analysis ranged from 1-5.

% example data used for an analysis so the code will work when tested
load data/Transects_LanduseClasses

for stream = startStream:endStream
    %clear out interior edge matrix and transectSamples each stream.
    interiorEdgeTemp=[];
    transectSamples_UniqueID_temp=[];

    %first clear out all temp matrices
    sampleXY_temp=[]; lidarGrdXY_temp=[]; lidarVegXY_temp=[]; lidarInteriorGrdXY_temp=[]; lidarInteriorVegXY_temp=[];

    
    %load lidar sample results for each mat file, this is used to plot final xy values
    dataName=(['results/lidarSampleResults_Strm',num2str(stream),'.mat']);
    load (dataName);
   
      
    % load sample data...
    dataName=(['results/forestSample_Strm',num2str(stream),'.mat']);
    load (dataName);
    dataTransect=(['results/transectWidth',num2str(stream),'.mat']);
    load (dataTransect);
    dataPoly=(['results/samplePolyVert',num2str(stream),'.mat']);
    load (dataPoly);
    dataName=(['results/lidarInteriorSampleResults_Strm',num2str(stream),'.mat']);
    load (dataName);
    dataPoly=(['results/sampleInteriorPolyVert',num2str(stream),'.mat']);
    load (dataPoly);
    %
    
    
    sampleXY_temp=(lidarSampleResults_temp); 
    forestSampleXY_temp=(forestSample);
   
    transectWidthXY_temp=(transectWidth_temp);
    samplePolyVertXY_temp=(samplePolyVert_StrmCurrent);
    sampleInteriorPolyVertXY_temp=(sampleInteriorPolyVert_StrmCurrent);
    
    if isempty(sampleXY)
        sampleXY=(sampleXY_temp); 
        
      %  forestSampleXY=(forestSampleXY_temp);
        TransectWidthXY=(transectWidthXY_temp);
        samplePolyVertXY=(samplePolyVert_StrmCurrent);
        LidarSampleResults=(lidarSampleResults_temp);
        LidarInteriorSampleResults=(lidarInteriorSampleResults_temp);
        sampleInteriorPolyVertXY=(sampleInteriorPolyVert_StrmCurrent);

    else
        sampleXY=[sampleXY; sampleXY_temp];

        TransectWidthXY=[TransectWidthXY; transectWidthXY_temp];
        samplePolyVertXY=[samplePolyVertXY; samplePolyVert_StrmCurrent];
        LidarSampleResults=[LidarSampleResults; lidarSampleResults_temp];
        LidarInteriorSampleResults=[LidarInteriorSampleResults; lidarInteriorSampleResults_temp];
        sampleInteriorPolyVertXY=[sampleInteriorPolyVertXY;sampleInteriorPolyVertXY_temp];

    end
    
    
%cycle through and assign edge, interior edge on their own row... also add
%unique id
m=0;

%first query the interior sample results for the current stream
LidarIntIndex_strmCur=find(LidarInteriorSampleResults(:,1)==stream);
LidarIntCurrent=LidarInteriorSampleResults(LidarIntIndex_strmCur,:);
%first query the interior sample results for the current stream
LidarIndex_strmCur=find(LidarSampleResults(:,1)==stream);
LidarCurrent=LidarSampleResults(LidarIndex_strmCur,:);

%then go through and add a unique ID column to the data
%%
%preallocated for transectSamples_
lenLidarCurrent=length(LidarCurrent(:,1));
transectSamples_UniqueID_temp=NaN(lenLidarCurrent,(length(LidarCurrent(1,:)) +1));

    for p=1:lenLidarCurrent
       %put together unique id string for each transect 
       strm=num2str(LidarCurrent(p,1));
       trans=num2str(LidarCurrent(p,2));
       direct=num2str(LidarCurrent(p,3));
       uniqueID=str2num(strcat(strm,trans,direct));
       transectSamples_UniqueID_temp(p,:)=[uniqueID LidarCurrent(p,:)];
       
    end 
    
    %1=interior, 2=bufferedge (landuse) 3= stream edge(stream)
    curSam=length(LidarIntCurrent(:,1));
    interiorEdgeTemp=NaN(curSam,17); %preallocate
    
    for p=1:(length(LidarIntCurrent(:,1)))
       %put together unique id string for each transect 
       strm=num2str(LidarIntCurrent(p,1));
       trans=num2str(LidarIntCurrent(p,2));
       direct=num2str(LidarIntCurrent(p,3));
       uniqueID=str2num(strcat(strm,trans,direct));
       %
       %add landuse class to interior / edge data
        %find the landuse code for current transect
        landuseType=[]; %clear out landuse Type
        landuseTypeIndex=(Transects_LanduseClasses(:,1)==uniqueID);
        landuseType=Transects_LanduseClasses(landuseTypeIndex,2);
        %if the transect isn't associated with a landuse, classify it as
        %999 (skip it in analysis)
        if isempty(landuseType)
            landuseType=999;
        end
       
       %}
           
       interiorEdgeTemp(p,:)=[uniqueID LidarIntCurrent(p,:) landuseType];
    end    
    
    %find transects with unique id's 
    uniqueTransectValues=unique(transectSamples_UniqueID_temp(:,1),'rows');
    
    % process each unique transect (each transect has multiple samples)
    
    for loop=1:length(uniqueTransectValues)
             
        %the current transect - 
        currentTransect=uniqueTransectValues(loop,1);
        %find all rows in the lidar samples matrix that have a value of
        %current transect
        allCurrentTransectRowsID=find(transectSamples_UniqueID_temp(:,1)==currentTransect);
        allCurrentTransectRows=transectSamples_UniqueID_temp(allCurrentTransectRowsID,:);
        
        %find the interior cell value for current transect
        interiorIndex=find((interiorEdgeTemp(:,1)==currentTransect) & (interiorEdgeTemp(:,8)== 1));
        landEdgeIndex=find((interiorEdgeTemp(:,1)==currentTransect) & (interiorEdgeTemp(:,8)== 2));
        
        %find the landuse code for current transect
        landuseType=[]; %clear out landuse Type
        landuseTypeIndex=(Transects_LanduseClasses(:,1)==currentTransect);
        landuseType=Transects_LanduseClasses(landuseTypeIndex,2);
        %if the transect isn't associated with a landuse, classify it as
        %999 (skip it in analysis)
        if isempty(landuseType)
            landuseType=999;
        end
        
        currentInter95th = (interiorEdgeTemp(interiorIndex,16));
        currentInter70th = interiorEdgeTemp(interiorIndex,15);  
        currentEdge70th =  interiorEdgeTemp(landEdgeIndex,15);  % not sure what this is or why i need it?? 
        currentInterFC = interiorEdgeTemp(interiorIndex,12);  
        currentEdgeFC =  interiorEdgeTemp(landEdgeIndex,12);
        
        % process all samples within each transect
        totalSamples=length(allCurrentTransectRows(:,1));
        max70th=max(allCurrentTransectRows(:,15));
        maxFC=max(allCurrentTransectRows(:,12));
        max10th=max(allCurrentTransectRows(:,13));
        max95th=(max(allCurrentTransectRows(:,16))); %this value was not converted when i processed the data so i am going to convert it to meters here.
        % add interior value to end cell 
        
        %preallocate
        finalTransects_temp=NaN(totalSamples,35);
        runningMeanFinal_Temp=NaN(totalSamples,8);
        for loop2=1:totalSamples
            
            % assign first value
            if loop2==1
                pct70Prev=0;
                pct70RelDiff=0;
                pct70CumRelDiff=0;
                runningMean=0;
               % pct70Current=0;
              %  70relDiffPrev=0;
            else         
            pct70Prev=allCurrentTransectRows((loop2-1),15);
            pct70Current=allCurrentTransectRows(loop2,15);
            
                if (loop2 ~= totalSamples) 
                pct70Next=allCurrentTransectRows(loop2+1,15);

                runningMean=((pct70Prev+pct70Current+pct70Next)/3);
                else
                    runningMean=0;
                end
            
                %calculate largest value
                if pct70Prev > pct70Current
                    largest70=pct70Prev;
                else
                    largest70=pct70Current;
                end
                
            % might add relative difference divided by the tallest veg
            pct70RelDiff=(pct70Current-pct70Prev) / largest70;          
            pct70CumRelDiff=pct70RelDiff+pct70CumRelDiff;
            %70relDiffPrev=pct70CumRelDiff;
            end
            
            if currentInter70th == 0
                diffValue70th =0;
                diffValueFC = 0;
            else
                diffValue70th = (currentInter70th - allCurrentTransectRows(loop2,15));
                diffValueFC = (currentInterFC - allCurrentTransectRows(loop2,12));
            end
            if loop2 ==totalSamples %deal with 999 value  
                streamEdgeLoc=999;%its the annoying fringe value that i calculated... 
            else
                streamEdgeLoc= (totalSamples-loop2-1)*5; %calculate location from stream edge
            end
            
            ratioLocation=(loop2/totalSamples);
            finalLocationRatio=(ratioLocation-((1/totalSamples)/2));
            
            % feel free to remove stats that aren't needed here
            finalTransects_temp(loop2,:) = [allCurrentTransectRows(loop2,:) currentInter70th currentEdge70th diffValue70th currentInterFC currentEdgeFC diffValueFC pct70RelDiff pct70CumRelDiff totalSamples loop2 ratioLocation finalLocationRatio maxFC max10th max70th streamEdgeLoc landuseType max95th currentInter95th]; %need to add max 95th to this
            
            runningMeanFinal_Temp(loop2,:)=[allCurrentTransectRows(loop2,1:5) allCurrentTransectRows(loop2,8) runningMean max70th];
        % transectSamples_UniqueID_temp
        end
        
        if isempty(FinalTransects)
            FinalTransects=finalTransects_temp; 
            RunningMeanFinal=runningMeanFinal_Temp;
        else
            FinalTransects=[FinalTransects; finalTransects_temp];
            RunningMeanFinal=[RunningMeanFinal; runningMeanFinal_Temp ];
        end
      
        
    end
    
    
    % now loop through each unique ID in the final results table and add
    % the interior value to the table. will then calculate a difference
    % value between the transect sample and the interior value.
    
   
    if isempty(InteriorEdge)      
        InteriorEdge=(interiorEdgeTemp);
    else
        InteriorEdge=[InteriorEdge; interiorEdgeTemp];
    end
 
end

%%
% reclassifiy landuse into two groups
 %landuse types 1:Ag; 2:Res 3:Infra 4:Urban Forest 5:Open Forest
landuseReclass=[1 1; 2 1 ; 3 1; 4 2; 5 2;999 999];
landuseReclass2=[1 1; 2 1 ; 3 2; 4 3; 5 4;999 999];

for k=1:length(InteriorEdge)
    newLanduse=[];
    newLanduse=landuseReclass(( landuseReclass(:,1) == InteriorEdge(k,17) ),2);
    newLanduse2=landuseReclass2(( landuseReclass2(:,1) == InteriorEdge(k,17) ),2);
    finalInteriorEdge(k,:)=[InteriorEdge(k,:) newLanduse newLanduse2]; %need to prellocate for this
end

% find unique transects
%calculate interior / ext/center sample point locations
    %1 = interior, 2=ME, landuse Edge (landuse) 3= stream edge(stream) 
    %4 = Matrix Fringe, 5 stream fringe,
    % type is column 8 in this matrix
    
    %fields 12-16 are the data
    % FC, 10th pct, 20thPct, 70thPt, 95thPct

  %remove transects where the landuse is 999  
  finalInteriorEdge_land = finalInteriorEdge((finalInteriorEdge(:,17) < 999),:);
   
 %remove transects where the width is too narrow (less than 23 m)
 finalInteriorEdge_Gr23 = finalInteriorEdge_land((finalInteriorEdge_land(:,5) > 24),:);
    
intEdgeUnique=unique(finalInteriorEdge_Gr23(:,1));

trans=length(intEdgeUnique(:,1));
for o=1:trans
    currentTrans=intEdgeUnique(o);
    %find inteEdgeTransects that have a unique id of the currentTrans
    currIntEdgeFinal=finalInteriorEdge_Gr23((finalInteriorEdge_Gr23(:,1)==currentTrans),:);
    currIntEdgeFinal2=finalInteriorEdge((finalInteriorEdge(:,1)==currentTrans),:);
    
    %calculate MEI for all data
    MEI_70_ME=funCalculateMEI(currIntEdgeFinal(2,15),currIntEdgeFinal(1,15));
    MEI_70_SE=funCalculateMEI(currIntEdgeFinal(3,15),currIntEdgeFinal(1,15));
    
    MEI_95_ME=funCalculateMEI(currIntEdgeFinal(2,16),currIntEdgeFinal(1,16));
    MEI_95_SE=funCalculateMEI(currIntEdgeFinal(3,16),currIntEdgeFinal(1,16));
    
    MEI_FC_ME=funCalculateMEI(currIntEdgeFinal(2,12),currIntEdgeFinal(1,12));
    MEI_FC_SE=funCalculateMEI(currIntEdgeFinal(3,12),currIntEdgeFinal(1,12));
    
    
    
    %fields: UniqueID, BufferWidth, 70thInt, 70thME, 70thSE,95thInt, 95thME,
    %95thSE,FCInt, FC_ME, FC_SE, 3 landuse classifications
    
    %this is used to create the final plots and determine how many points to
    %plot.
    widthClasses=[23 30 2330; 31 50 3050; 51 100 50100; 101 150 100149; 151 199 150199; 200 200 200];
    % loop through to figure out what width class each transect is in
    widthClassFinal=NaN;
    for l=1:length(widthClasses(:,1))
        
        widthStart=widthClasses(l,1);
        widthEnd=widthClasses(l,2);
        
        if (currIntEdgeFinal(1,5) >= widthStart) && (currIntEdgeFinal(1,5)<= widthEnd)
            widthClassFinal=widthClasses(l,3); 
        end
    end
    % fields as follows
    % UniqueID, Width, WidthClass, 70thPctInterior, 70thME,
    % 70thSE,95thPctInterior, 95thME, 95thSE, FCInterior, FCME, FC_SE
    % MEI_70_ME,MEI_70_SE,MEI_95_ME,MEI_95_SE,MEI_FC_ME,MEI_FC_SE, 
    % landuse 1, landuse 2, landuse 3,
    
    EdgeIntOutput(o,:)=[currIntEdgeFinal(1,1),currIntEdgeFinal(1,5),widthClassFinal,currIntEdgeFinal(1,15),currIntEdgeFinal(2,15),currIntEdgeFinal(3,15), currIntEdgeFinal(1,16),currIntEdgeFinal(2,16),currIntEdgeFinal(3,16), currIntEdgeFinal(1,12),currIntEdgeFinal(2,12),currIntEdgeFinal(3,12),MEI_70_ME,MEI_70_SE,MEI_95_ME,MEI_95_SE,MEI_FC_ME,MEI_FC_SE,currIntEdgeFinal(1,17:19)];
    EdgeInt_MEIOut(o,:)=[currIntEdgeFinal(1,1),currIntEdgeFinal(1,5),widthClassFinal, MEI_70_ME,MEI_70_SE,MEI_95_ME,MEI_95_SE,MEI_FC_ME,MEI_FC_SE];
end


%%

%create a final interior / edge matrix that i will use for analysis

% cycle through each interior / edge row and pull out the following fields
% unique ID, 70thPct_int,70thPct_ME,70thPct_SE, 95thPct_int,95thPct_ME,
% 95thPct_SE, FC_Int,FC_ME,FC_SE, Landuse, landuse 2, landuse 3



% SAVE RESULTS
% ??, totalGroundPoints, totalVegPoints, VegPts greater than 2m, Fractional
% Cover, 10, 20, 70, 90 (percentiles)
%%
%save out final transects as it will be used to calculate MEI values
save results/finalTransects/FinalTransects FinalTransects

save results/finalTransects/InteriorEdge InteriorEdge

% export InteriorEdge at text file
%%
%calculate interior / ext/center sample point locations
    %1 = interior, 2=bufferedge (landuse) 3= stream edge(stream) 
    %4 = Matrix Fringe, 5 stream fringe,
    
    
    %UniqueID, Width, WidthClass, 70thPctInterior, 70thME,
    % 70thSE,95thPctInterior, 95thME, 95thSE, FCInterior, FCME, FC_SE
    % MEI_70_ME,MEI_70_SE,MEI_95_ME,MEI_95_SE,MEI_FC_ME,MEI_FC_SE, 

    fnam=('results/finalTransects/InterEdge_Final.txt');
     hdr={'Unique ID,','Width,','WidthClass,','70pctInt,','70pctME,','70pctSE,','95pctInt,','95pctME,','95pctSE,','FC_Int,','FC_ME,','FC_SE,','MEI_70_ME,','MEI_70_SE,','MEI_95_ME,','MEI_95_SE,','MEI_FC_ME,','MEI_FC_SE,','landuse1,','landuse2,','landuse3,'};
     txt=sprintf('%s',hdr{:});
     txt(end)='';
     dlmwrite(fnam,txt,'');
     dlmwrite(fnam,EdgeIntOutput,'-append','precision','%8.2f','delimiter',',');
     
    %{
     fnam=('results/all/InteriorEdgeResults.txt');
     hdr={'Unique ID,','StreamName,','TransectNum,','pos/neg,','BufWidth,','xCoord,','YCoord,','EdgeInterPos,','totalGrdPts,','totalVegPts,','VegPtsGr2m,','FC,','10th,','20th,','70th,','95th,','landuse,'};
     txt=sprintf('%s',hdr{:});
     txt(end)='';
     dlmwrite(fnam,txt,'');
     dlmwrite(fnam,InteriorEdge,'-append','precision','%8.2f','delimiter',',');
%}     

     fnam=('results/finalTransects/allFinalTransects.txt');
     hdr={'UniqueID,','StreamName,','TransectNum,','Pos_Neg,','BufWidth,','xCoord,','YCoord,','DistFromEdge,','totalGrdPts,','totalVegPts,','VegPtsGr2m,','FC,','10th,','20th,','70th,','95th,','70th_Interior,','70thLandEdge,','Diff70_IntMinusEdge,','FCInt,','FCEdge,','FCDiff_InterMinusEdge,','70threlDiff,','70thCumRelDiff,','totalSamples,','sampleNumber,','RatioLocation,','FinalLocationRatio,','MaxFC,','Max10,','Max70,','DistFromStream Edge,','Landuse Class,','max95ThPct,','95thInt,'};
     txt=sprintf('%s',hdr{:});
     txt(end)='';
     dlmwrite(fnam,txt,'');
     dlmwrite(fnam,FinalTransects,'-append','precision','%8.2f','delimiter',',');
     

% write out the initial transect results (with the buffer width, and fractional cover for the entire transect (5m x xwidth) with headers
     fnam=('results/finalTransects/transectResults.txt'); % <-  data file
     hdr={'OrigStrm,','StreamName,','TransectNum,','Pos/Neg,','X_Stream,','Y_Stream,','X_end,','Y_end,','Width,','pctCover,','X_begin,','Y_begin,','StreamWidth'};
     txt=sprintf('%s',hdr{:});
     txt(end)='';
     dlmwrite(fnam,txt,'');
     dlmwrite(fnam,TransectWidthXY,'-append','precision','%8.2f','delimiter',',');
    
   
%Uncomment the code below to plot the final transects to make sure they look correct.
     %{
     
% Plot the interior data...
figure(1)

%
plot (sampleInteriorPolyVertXY(:,14), sampleInteriorPolyVertXY(:,15), 'ko', 'MarkerFaceColor','k','MarkerSize',3), hold on
plot (sampleInteriorPolyVertXY(:,8), sampleInteriorPolyVertXY(:,9), 'ko', 'MarkerFaceColor','k','MarkerSize',3), hold on
plot (sampleInteriorPolyVertXY(:,10), sampleInteriorPolyVertXY(:,11), 'ko', 'MarkerFaceColor','k','MarkerSize',3), hold on
plot (sampleInteriorPolyVertXY(:,12), sampleInteriorPolyVertXY(:,13), 'ko', 'MarkerFaceColor','k','MarkerSize',3), hold on


% plot the sample location (center point) of each sample box
plot(sampleXY(:,5), sampleXY(:,6),'rx', 'MarkerSize',4), hold on

plot (lidarInteriorGrdXY(:,1),lidarInteriorGrdXY(:,2),'k*'), hold on
plot (lidarInteriorVegXY(:,1),lidarInteriorVegXY(:,2),'go', 'MarkerFaceColor','g', 'MarkerSize',3), hold on

     
for f=startStream:endStream
    vertices=AllStreamLength(f,3);
    A2=zeros(length(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
axis equal

figure(2)

%plot (x1,y1,'gs', x2,y2,'rs', x3,y3,'bs', x4,y4,'ks')
%Plot the Lidar Data used to calculate the xm by xm sample box values

plot (lidarGrdXY(:,1),lidarGrdXY(:,2),'k*'), hold on
plot (lidarVegXY(:,1),lidarVegXY(:,2),'go', 'MarkerFaceColor','g', 'MarkerSize',3), hold on

% plot the sample location (center point) of each sample box
plot(sampleXY(:,5), sampleXY(:,6),'ro', 'MarkerSize',10), hold on

% plot beginning and end of buffer transect (based upon final calculated width
plot(TransectWidthXY(:,5), TransectWidthXY(:,6),'go', 'MarkerFaceColor','g','MarkerSize',10), hold on
plot(TransectWidthXY(:,7), TransectWidthXY(:,8),'m*', 'MarkerSize',10), hold on
plot(TransectWidthXY(:,11), TransectWidthXY(:,12),'m*', 'MarkerSize',10), hold on

%
plot (samplePolyVertXY(:,14), samplePolyVertXY(:,15), 'ko', 'MarkerFaceColor','k','MarkerSize',3), hold on
plot (samplePolyVertXY(:,8), samplePolyVertXY(:,9), 'ko', 'MarkerFaceColor','k','MarkerSize',3), hold on
plot (samplePolyVertXY(:,10), samplePolyVertXY(:,11), 'ko', 'MarkerFaceColor','k','MarkerSize',3), hold on
plot (samplePolyVertXY(:,12), samplePolyVertXY(:,13), 'ko', 'MarkerFaceColor','k','MarkerSize',3), hold on

%plot (LidarSampleResults(:,5),LidarSampleResults(:,6), 'kx'), hold on
allStreamsLen=length(AllStreamLength);

%for f=1:allStreamsLen
 for f=startStream:endStream
    vertices=AllStreamLength(f,3);
    A2=zeros(length(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
axis equal

%}

disp('all done-- all data exported!!');