/
Mapping aboveground biomass from individual tree segmentation data

Mapping aboveground biomass from individual tree segmentation data

Jan 09, 2016

Variable Area Local Maxima

I generated stem maps for each LiDAR tile from the Boulder Creek Snow-off dataset using my local maxima stem isolation algorithm (Swetnam and Falk 2014). Data are available from OpenTopography.org

The algorithm has been updated to include an inverse watershed segmentation (right panel) which reports the minor and major axis of canopy diameter and the total canopy area. The equivalent diameter of the canopy (min+max/2) is shown in the left panel (yellow circles). Some watersheds are empty if there was no tree identified within it (this can be caused by deformed canopy shapes that result in the watershed segmentation splitting the canopy up in to multiple areas. Conversely, some watersheds may have more than one tree identified within them, e.g. a very large tree and a smaller seedling tree off the edge of its canopy. 

The code for the VLM:

% Variable area local maxima (VLM) search method Version 1.1,  
% Date: 7/14/2015 
% Author: Tyson L. Swetnam, University of Arizona. Email:
% tswetnam@email.arizona.edu,
% Version 1.0 was published in 
% Swetnam, T. L., and Falk, D. A. (2014)  
% Application of Metabolic Scaling Theory to reduce error in local maxima
% tree segmentation from aerial LiDAR. Forest Ecology and Management, 323,
% 158-167. http://dx.doi.org/10.1016/j.foreco.2014.03.016,
%
%
% This algorithm outputs an individual tree identifier number (ID), the X
% and Y coordinates in units of pixels - these are reprojected in spatial,
% units using the export_utm.m option, the tree height (Z),
% and predicted canopy radius (PRED), canopy area in units of the data
% (AREA), canopy equivalent diameter (EQDIAM), major axis length
% (MAJAXLGTH), minor axis length (MINAXLGTH), maximum height of inverse
% watershed (MAXHT), minimum height of inverse watershed (MINHT), and mean
% hieght of pixels in the inverse watershed (MEANHT)
% Data are set in units of the projected ASC or TIF.
% The intended purpose for this algorithm is locating local maxima (trees),  
% derived from an aerial LiDAR Canopy Height Model.
% This script first requires an input Matlab structure file generated by
% the asc_import.m. The VLM first conducts an anisotropic diffusion of the
% surface and dilates the image to reduce no-data pixels. Next, it isolates
% the local maxima and collates them into a single file index. The second
% part of the execution involves a k-nearest neighbor exhaustive searcher
% which uses a logic rule to eliminate individual stems within the minimum
% canopy radius of a larger tree. This local maxima code was adapted from
% work found on the MATLAB Central File Exchange published by Yonathan
% Nativ:
% http://www.mathworks.com/matlabcentral/fileexchange/14498-local-maxima-minima,
% Anisodiff2D.m was written by Daniel Lopes and is available on the 
% MATLAB File Exchange 
% http://www.mathworks.com/matlabcentral/fileexchange/14995-anisotropic-diffusion-perona-malik/content/anisodiff_Perona-Malik/anisodiff2D.m
% the purpose of the Anisotropic diffusion is to Gaussian smooth data
% gaps in the canopy height model wile retaining edge definition.
function stems = vlm(img_import,canopyratio,anisotropiciterations)
% Converts a uint8 (e.g. TIF,JPG) file into a double, note that TIF do not
% always have data in height units, rather data can be coded as 8-bit 0-255
% versus decimal floating point values
tic
img=img_import.noneg;
[r1, c1]=size(img);
% Add the border thicknesses to the total image size
r2 = r1 + 2.*100;
c2 = c1 + 2.*100;
img2=repmat(eps, [r2 c2 1]);
img2((100+1):(r2-100), (100+1):(c2-100),:) = img;
x1=double(img2);
% Creates grid 'xold', this retains the grid's original values.
xold=x1;
% 2D anisotropic diffusion with 's' iterations blurs the surface while
% retaining edge values
% the greater the number of iterations the smoother the output surface, 
% the less likely smaller maxima are retained near larger ones
% Settings for the anisodiff2D are set to default settings from original
% (anisodiff2d.m):
% x1 = grayscale image or floating point values, s = number of iterations, 
% maximum Delta_T = 1/7, Kappa = 30, 
% Option = 2 privileges wide regions over small.
x2new=anisodiff2D(x1,anisotropiciterations,1./7,30,2);
SE=strel('disk',1,0);
% Creates an elliptical shaped neighborhood size that functions as
% the limiting area of the variable window for any other local maxima to be
% recorded, forces favoring of larger trees over smaller ones as disk size
% increases, 
X2=imdilate(x2new,SE);
% Output becomes the X and Y coordinates for row# and column# of stems >h 
[r2, c2]=find(x2new==X2&xold>2);
% Output to determine stem heights
f2=find(x2new==X2&xold>2);
% Trees must be greater than 2 units (intended to be in meters)
% Final output file, should be 3 columns = X,Y,Z
stems2=[c2,r2,x1(f2)];
stemsall=cat(1,stems2);
x=sortrows(stemsall,3);
[u,m,n] = unique(x,'rows'); % Determine unique values 
A=accumarray(n(:),1); % Finds any duplicate points and removes them
stem=cat(2,u,A);
% Outputs to structure
stem1=stem;
stems.x=stem1(:,1); % X-axis in pixel number to be converted to UTME
stems.y=stem1(:,2); % Y-axis in pixel number to be converted to UTMN
stems.z=stem1(:,3); % stem height in units of the data
stems.id=transpose(1:(nnz(stems.x)));
% Self-thins the number of stems by their inter-stem distances 
% Nearest neighbor searcher attempts to eliminate errors of commission
% based on expected minimum canopy radius from stems.r
NS1=ExhaustiveSearcher([stems.x stems.y]);
stems.pred=canopyratio.*stems.z; 
% normalization constant for canopy radius called by the self-thinning
% algorithm, 
[idx1,dist1]=rangesearch(NS1,[stems.x stems.y],(8./img_import.cellsize)); 
% set to spherical radius of the tree
list_id_flagged1=zeros(10,size(stems.y,1));
for i=1:1:size(stems.x,1); 
    stems.z(i)=stems.z(i)+(rand(1,1)./1000); 
    % adds a tiny random number to each height to eliminate the chance of
    % two equal height maxima 
end
for i=1:1:size(stems.x,1);
   	neighbor_ID1=transpose(cell2mat(idx1(i,1)));
	neighbor_dist1=transpose(cell2mat(dist1(i)));
    neighbor_z=stems.z(neighbor_ID1);
    flag1= neighbor_dist1 < (stems.pred(i)./img_import.cellsize) & neighbor_z < stems.z(i);
    list_id_flagged1(1:size(neighbor_ID1(flag1),1),i)=neighbor_ID1(flag1);
end
C1=list_id_flagged1;
B1=unique(C1);
B1(1,:)=[]; % deletes first sorted value [a zero] values are ID
A1.id=B1(:,1); %ID
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Recompiles the tree list, removing commission error stems
stems.lia=ismember(stems.id,A1.id);
for j=1:1:size(stems.x,1);
    if stems.lia(j,1)==0;
        stems.x(j,1)=stems.x(j,1);
        stems.y(j,1)=stems.y(j,1);
        stems.z(j,1)=stems.z(j,1);
        stems.pred(j,1)=stems.pred(j,1);  
    elseif stems.lia(j,1)==1;
        stems.x(j,1)=NaN;
        stems.y(j,1)=NaN;
        stems.z(j,1)=NaN;
        stems.pred(j,1)=NaN;
    end
   
end
STEMS = [stems.x stems.y stems.z stems.pred];
STEMS = STEMS(~any(isnan(STEMS),2),:);
stems.x = STEMS(:,1)-100;
stems.y = STEMS(:,2)-100;
stems.z = STEMS(:,3);
stems.pred = STEMS(:,4);
stems.id=transpose(1:(nnz(stems.x)));
% Inverse watershed segmentation for canopy radius and area estimate
for i=1:size(stems.x);
stems.watershed_id(i,1)=img_import.watershed(stems.y(i),stems.x(i));
end
stems.STATS=regionprops(img_import.watershed,img_import.noneg,'Area','EquivDiameter','MajorAxisLength','MinorAxisLength','MaxIntensity','MinIntensity','MeanIntensity');
stems.watershed_index=double(stems.watershed_id);
for i=1:size(stems.x);
    if stems.watershed_index(i)>0;
        stems.area(i,1)=stems.STATS(stems.watershed_index(i)).Area;
        stems.equivdiam(i,1)=stems.STATS(stems.watershed_index(i)).EquivDiameter;
        stems.majaxlgth(i,1)=stems.STATS(stems.watershed_index(i)).MajorAxisLength;
        stems.minaxlgth(i,1)=stems.STATS(stems.watershed_index(i)).MinorAxisLength;
        stems.maxht(i,1)=stems.STATS(stems.watershed_index(i)).MaxIntensity;
        stems.minht(i,1)=stems.STATS(stems.watershed_index(i)).MinIntensity;
        stems.meanht(i,1)=stems.STATS(stems.watershed_index(i)).MeanIntensity;
    else
        stems.area(i,1)=NaN;
        stems.equivdiam(i,1)=NaN;
        stems.majaxlgth(i,1)=NaN;
        stems.minaxlgth(i,1)=NaN;
        stems.maxht(i,1)=NaN;
        stems.minht(i,1)=NaN;
        stems.meanht(i,1)=NaN;
    end
end
stems.area=stems.area.*(img_import.cellsize.^2);
stems.equivdiam=stems.equivdiam.*img_import.cellsize;
stems.majaxlgth=stems.majaxlgth.*img_import.cellsize;
stems.minaxlgth=stems.majaxlgth.*img_import.cellsize;
toc

After I generated the stem maps for each tile I merged the tiles, reordered the tree list by UTM easting, and assigned new identifier numbers:

tile=asc_import('I:\BC2010_snoff\Catchment_GIS\Betasso\Canopy_Height\ot_470500_4428000.asc');
tile_v=vlm(tile,0.15,2);
tile_exp=export_utm(tile_v,tile);
dlmwrite('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_470500_4428000.csv',tile_exp,'delimiter',',','precision',10);
tile=asc_import('I:\BC2010_snoff\Catchment_GIS\Betasso\Canopy_Height\ot_470500_4429000.asc');
tile_v=vlm(tile,0.15,2);
tile_exp=export_utm(tile_v,tile);
dlmwrite('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_470500_4429000.csv',tile_exp,'delimiter',',','precision',10);
tile=asc_import('I:\BC2010_snoff\Catchment_GIS\Betasso\Canopy_Height\ot_471000_4428000.asc');
tile_v=vlm(tile,0.15,2);
tile_exp=export_utm(tile_v,tile);
dlmwrite('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_471000_4428000.csv',tile_exp,'delimiter',',','precision',10);
tile=asc_import('I:\BC2010_snoff\Catchment_GIS\Betasso\Canopy_Height\ot_471000_4429000.asc');
tile_v=vlm(tile,0.15,2);
tile_exp=export_utm(tile_v,tile);
dlmwrite('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_471000_4429000.csv',tile_exp,'delimiter',',','precision',10);
tile=asc_import('I:\BC2010_snoff\Catchment_GIS\Betasso\Canopy_Height\ot_471500_4428000.asc');
tile_v=vlm(tile,0.15,2);
tile_exp=export_utm(tile_v,tile);
dlmwrite('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_471500_4428000.csv',tile_exp,'delimiter',',','precision',10);
tile=asc_import('I:\BC2010_snoff\Catchment_GIS\Betasso\Canopy_Height\ot_471500_4429000.asc');
tile_v=vlm(tile,0.15,2);
tile_exp=export_utm(tile_v,tile);
dlmwrite('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_471500_4429000.csv',tile_exp,'delimiter',',','precision',10);
tic
ot_470500_4428000=dlmread('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_470500_4428000.csv',',');
ot_470500_4429000=dlmread('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_470500_4429000.csv',',');
ot_471000_4428000=dlmread('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_471000_4428000.csv',',');
ot_471000_4429000=dlmread('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_471000_4429000.csv',',');
ot_471500_4428000=dlmread('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_471500_4428000.csv',',');
ot_471500_4429000=dlmread('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\ot_471500_4429000.csv',',');
merged=vertcat(ot_470500_4428000,ot_470500_4429000,ot_471000_4428000,ot_471000_4429000,ot_471500_4428000,ot_471500_4429000);
merged=sortrows(merged,2);
n=numel(merged);
myvector=[1:n];
dlmwrite('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\Betasso_stems.csv',merged,'delimiter',',','precision',10);
col_header={'ID','UTME','UTMN','HT','PRED','AREA','EQDIAM','MAJAX','MINAX','MAXHT','MINHT','MEANHT'};
xlswrite('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\Betasso_stems.xlsx',merged,1,'A2');
xlswrite('I:\BC2010_snoff\Catchment_GIS\Betasso\Stem_Map\Betasso_stems.xlsx',col_header,1,'A1');
toc

After I created the tiled data, I clipped the stems to each of the three catchments based on the GIS shape file polygon.

Estimating Individual tree aboveground biomass and carbon

To calculate the biomass of each individual tree I searched the literature for allometric models based on tree height and canopy diameter. There is not much research to date on tree height-to-biomass calculation alone, so I am using my published diameter at breast height (DBH) prediction model from Swetnam and Falk (2014):

DBH (cm) = Beta * HT (m) * √CanopyDiameter (m); RMSE = 7.66 cm, r2 = 0.811

I can estimate the above ground biomass/Carbon (AGB/AGC) using published wood specific gravity (SG, 1Chojnacky et al. 2014) and C content (C%, 2Lamlom and Savidge 2003). Biomass equations are from Chojnacky et al. (2014) where the form of the equation is: ln(AGB) = β0 + β1ln(DBH).

SpeciesSpecific Gravity (SG)1Carbon Density (%)2β0 β1R2-statistic1
Abies concolor0.37 -3.17742.64260.75
Abies lasiocarpa0.31 -2.31232.34820.75
Abies spp. 48.55 - 50.08   
Picea engelmannii0.33 -3.03002.55670.81
Picea spp.0.3749.95 - 50.39   
Pinus ponderosa0.3852.47 ± 0.38-2.61772.4638

0.83

Pinus flexilis0.36    
Pinus contorta0.3850.32 ± 0.43-2.61772.46380.83
Populus tremuloides0.3547.09 ± 0.75   
Pseudotsuga menziesii 0.4550.50 ± 0.36-2.46232.48520.86
Quercus gambellii0.61    

Lodgepole pine (Pinus contorta) allometry is reported by Litton et al. (2004) as:  (height = 2.89 * DBH^0.59, r2 = 0.82, P= 0.01, n = 200). The inverse model for height based DBH is: 0.1655 * Ht ^ 1.695. The derived model for AGB based on height assuming the pipe-model, and a wood specific gravity SG = 0.38 for lodgepole pine from Chojnacky et al. (2014): AGB (kg) = 0.0008175 * Ht (m) ^ 4.39 

As an example of the estimated above ground carbon (AGC) derived from the Chojnacky et al. (2014) Pinus spp. DBH based AGB model for P. ponderosa and P. contorta (same model) versus the pipe model which uses height and DBH (as radius), I graph the Betasso catchment trees below, note that both the X and Y axis use the same estimated DBH values. 

I used the LANDFIRE Existing Vegetation Type (EVT) to characterize individual tree forest types. By catchment the dominant EVTs are:

Betasso

  • Southern Rocky Mountain Ponderosa Pine Woodland (98.166%)
  • Developed or Ruderal Shrubland (0.765%)
  • All Other Classes (1.1%)

Gordon Gulch

  • Southern Rocky Mountain Ponderosa Pine Woodland (41.923%)
  • Rocky Mountain Lodgepole Pine Forest (41.669%)
  • Inter-Mountain Basins Aspen-Mixed Conifer Forest and Woodland (4.338%)
  • Southern Rocky Mountain Dry-Mesic Montane Mixed Conifer Forest and Woodland (4.064)
  • Rocky Mountain Montane Riparian Forest and Woodland (2.933%)
  • Southern Rocky Mountain Mesic Montane Mixed Conifer Forest and Woodland (2.660%)
  • Rocky Mountain Aspen Forest and Woodland (2.042%)
  • All Other Classes (0.371%)

Como Creek

  • Rocky Mountain Subalpine Dry-Mesic Spruce-Fir Forest and Woodland (57.384%)
  • Rocky Mountain Lodgepole Pine Forest (17.384%)
  • Rocky Mountain Subalpine Mesic-Wet Spruce-Fir Forest and Woodland (16.965%)
  • Inter-Mountain Basins Aspen-Mixed Conifer Forest and Woodland (2.429%)
  • Snow-Ice (4.084%)
  • All Other Classes (1.754%)

Table: The EVT's by Catchment and above ground biomass/carbon models available for each.

LANDFIRE EVT CodeVegetation Class Descriptions

Betasso

(n, 30 m2)

Betasso

(% of Area)

Gordon Gulch

(n, 30 m2)

Gordon Gulch

(% of Area)

Como Creek

(n, 30 m2)

Como Creek

(% of Area)

Height-to-Biomass Allometry Model ParametersRMSE (kg)R2Reference
3011Rocky Mountain Aspen Forest and Woodland10 0.05836462.04213490.425AGB (kg) = 0.006 * ht (m) ^ 3.947643.50.492Swetnam 2013
3049Rocky Mountain Foothill Limber Pine-Juniper Woodland  180.010      
3050Rocky Mountain Lodgepole Pine Forest  7436541.6695515417.384   Lefsky et al. 1999
3051Southern Rocky Mountain Dry-Mesic Montane Mixed Conifer Forest and Woodland  72534.064  AGB (kg) = 0.027 * ht (m) ^ 3.605668.20.717Swetnam 2013
3052Southern Rocky Mountain Mesic Montane Mixed Conifer Forest and Woodland  47482.660  AGB (kg) = 0.039 * ht (m) ^ 3.455779.80.755Swetnam 2013
3054Southern Rocky Mountain Ponderosa Pine Woodland1696898.1667481941.923  AGB (kg) = 0.010 * ht (m) ^ 3.935526.70.760Swetnam 2013
3055 Rocky Mountain Subalpine Dry-Mesic Spruce-Fir Forest and Woodland  790.044218182757.384    
3056Rocky Mountain Subalpine Mesic-Wet Spruce-Fir Forest and Woodland    5382516.965    
3057Rocky Mountain Subalpine-Montane Limber-Bristlecone Pine Woodland   1080.0611710.054    
3061Inter-Mountain Basins Aspen-Mixed Conifer Forest and Woodland10 0.05877424.33877062.429AGB (kg) = 0.027 * ht (m) ^ 3.605668.20.717Swetnam 2013
3070Rocky Mountain Alpine Dwarf-Shrubland     830.026    
3080Inter-Mountain Basins Big Sagebrush Shrubland  1230.069      
3086Rocky Mountain Lower Montane-Foothill Shrubland  2200.123400.013    
3092Southern California Coastal Scrub    400.013    
3144Rocky Mountain Alpine Turf    9160.289N/A   
3145Rocky Mountain Subalpine-Montane Mesic Meadow    270.009N/A   
3159Rocky Mountain Montane Riparian Forest and Woodland  52352.933      
3182Introduced Upland Vegetation-Perennial Grassland and Forbland  610.034  N/A   
3219Inter-Mountain Basins Sparsely Vegetated Systems II    1080.034    
3220Artemisia tridentata ssp. vaseyana Shrubland Alliance180.104  190.006    
3222Rocky Mountain Alpine/Montane Sparsely Vegetated Systems II     240.008    
3252 Rocky Mountain Subalpine/Upper Montane Riparian Shrubland   510.02930250.953    
3292Open Water      --   
3293Snow-Ice    129574.084--   
3297Developed-Medium Intensity1320.764    --   
3900Western Cool Temperate Urban Deciduous Forest 130.075        
3923Western Cool Temperate Developed Ruderal Shrubland 340.197        
3924Western Cool Temperate Developed Ruderal Grassland    20.0006N/A   
SG=[];
C=[];
for i=1:size(EVT);
    if EVT(i) ==3011 % Rocky Mountain Aspen Forest and Woodland
        SG(i,1) = 0.35;
        C(i,1) = 0.47;
    elseif EVT(i) == 3049 % Rocky Mountain Foothill Limber Pine-Juniper Woodland
        SG(i,1) = 0.38;
        C(i,1) = 0.5;
    elseif EVT(i) == 3050 % Rocky Mountain Lodgepole Pine Forest
        SG(i,1) = 0.38;
        C(i,1) = 0.5;
    elseif EVT(i) == 3051 % Southern Rocky Mountain Dry-Mesic Montane Mixed Conifer Forest and Woodland
        SG(i,1) = 0.38;
        C(i,1) = 0.5;
    elseif EVT(i) == 3052 % Southern Rocky Mountain Mesic Montane Mixed Conifer Forest and Woodland
        SG(i,1) = 0.38;
        C(i,1) = 0.5;
    elseif EVT(i) == 3054 % Southern Rocky Mountain Ponderosa Pine Woodland
        SG(i,1) = 0.38;
        C(i,1) = 0.52;
    elseif EVT(i) == 3055 % Rocky Mountain Subalpine Dry-Mesic Spruce-Fir Forest and Woodland
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 3056 % Rocky Mountain Subalpine Mesic-Wet Spruce-Fir Forest and Woodland
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 3057 % Rocky Mountain Subalpine-Montane Limber-Bristlecone Pine Woodland 
        SG(i,1) = 0.38;
        C(i,1) = 0.5;
    elseif EVT(i) == 3061 % Inter-Mountain Basins Aspen-Mixed Conifer Forest and Woodland
        SG(i,1) = 0.38;
        C(i,1) = 0.5;
    elseif EVT(i) == 3070 % Rocky Mountain Alpine Dwarf-Shrubland 
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 3080 % Inter-Mountain Basins Big Sagebrush Shrubland
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 3086 % Rocky Mountain Lower Montane-Foothill Shrubland
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 3092 % Southern California Coastal Scrub
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 3144 % Rocky Mountain Alpine Turf
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 3145 % Rocky Mountain Subalpine-Montane Mesic Meadow
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 3159 % Rocky Mountain Montane Riparian Forest and Woodland
        SG(i,1) = 0.33;
        C(i,1) = 0.5;  
    else %Non Forested - missclassified
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    end
end

 

Table: Muldavin et al. (2006) Vegetation Map Types

VALUECOUNTCLASS_NAMEHectares (ha)SGCLa JaraLa JaraJaramilloJaramilloHistory GroveHistory Grove  
0 NaN 0.330.5        
14353949Spruce-Fir Forest and Woodland (Dry Mesic)10633.970.330.5        
22732982Spruce-Fir Forest and Woodland (Moist Mesic)6674.9710.330.5        
32320514Forest Meadow5667.5770.330.5        
422084678Mixed Conifer Forest and Woodland (Dry Mesic)53939.10.380.5        
514126495Mixed Conifer Forest and Woodland (Moist Mesic)34502.140.380.5        
7783518Blue Spruce Fringe Forest1913.6480.330.5        
103241666Aspen Forest and Woodland (Dry Mesic)7917.360.350.47        
111921226Aspen Forest and Woodland (Moist Mesic)4692.3570.350.47        
139348740Ponderosa Pine Forest22833.130.380.52        
141459721Gambel Oak-Mixed Montane Shrubland3565.1920.610.5        
164990766Upper Montane Grassland12189.320.330.5        
1712778582Lower Montane Grassland31209.980.330.5        
195900099Wet Meadow14410.270.330.5        
201033140Wetland2523.3270.330.5        
2114647Montane Riparian Shrubland35.780660.330.5        
24161175Sparsely Vegetated Rock Outcrop393.66140.330.5        
25925811Felsenmeer Rock Field2261.1750.330.5        
261554344Roads-Disturbed Ground3796.2830.330.5        
2756340Open Water137.61220.330.5        
2816769Post-Fire Bare Ground40.945130.330.5        

 

Matlab Script for EVT classification of Jemez_stems_evt.csv imported into Matlab

SG=[size(EVT),1];
C=[size(EVT),1];
for i=1:size(EVT);
    if EVT(i) <=3 %NaNs, Spruce-Fir Forest and Woodland (Dry Mesic), Spruce-Fir Forest and Woodland (Moist Mesic), & Forest Meadow
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 4 %Mixed Conifer Forest and Woodland (Dry Mesic)
        SG(i,1) = 0.38;
        C(i,1) = 0.5;
    elseif EVT(i) == 5 %Mixed Conifer Forest and Woodland (Moist Mesic)	
        SG(i,1) = 0.38;
        C(i,1) = 0.5;
    elseif EVT(i) == 7 %Blue Spruce Fringe Forest
        SG(i,1) = 0.33;
        C(i,1) = 0.5;
    elseif EVT(i) == 10 %Aspen Forest and Woodland (Dry Mesic)
        SG(i,1) = 0.33;
        C(i,1) = 0.47;
    elseif EVT(i) == 11 %Aspen Forest and Woodland (Moist Mesic)
        SG(i,1) = 0.33;
        C(i,1) = 0.47;
    elseif EVT(i) == 13 %Ponderosa Pine Forest
        SG(i,1) = 0.38;
        C(i,1) = 0.52;
    elseif EVT(i) == 14 %Gambel Oak-Mixed Montane Shrubland
        SG(i,1) = 0.61;
        C(i,1) = 0.5;
    else %Non Forested - missclassified
        SG(i,1) = 0.33;
    end
end

Next I estimate the DBH of each tree based on the tree height multiplied by the square root of its estimated canopy diameter (as measured by the VLM's inverse watershed segmentation). If the tree had an estimated canopy diameter that was >150% of its predicted canopy diameter the predicted canopy diameter is used instead. This is because the watershed segmentation does not always represent a canopy properly, in particular for smaller trees next to very large trees. The normalization constant is estimated from data to be ~0.82 for all trees (data are from a mix of Pseudotsuga, Picea, and Pinus trees measured in Arizona and New Mexico, Swetnam and Falk 2014). 

% If the equivalent diameter (major axis+minor axis/2) is >300% of the predicted minimum canopy radius
% the diameter defaults to the predicted canopy diameter. 
CD=1:size(EVT);
for i =1:size(EQDIAM);
    if EQDIAM(i) > (3.*PRED(i));
        CD(i)=PRED(i);
    else
        CD(i)=EQDIAM(i);
    end
end
DBH=1:size(EVT);
AGB=1:size(EVT);
AGC=1:size(EVT);
%Assumes that HT and CD are in units of meters, DBH is in units centimeters
for i=1:size(EVT);
    DBH(i)=0.82.*(HT(i).*sqrt(CD(i)));
    AGB(i)=(HT(i).*100).*3.1415952.*((DBH(i)./2).^2).*SG(i)./1000; %Calculate aboveground biomass (AGC) in units of Kg
    AGC(i)=AGB(i).*C(i);
end
DBH=transpose(DBH);
AGC=transpose(AGC);
AGB=transpose(AGB);
CD=transpose(CD);

For an EVT specific classification I tried:

% If the equivalent diameter (major axis+minor axis/2) is >300% of the predicted minimum canopy radius
% the diameter defaults to the predicted canopy diameter. 
CD=1:size(EVT);
for i =1:size(EQDIAM);
    if EQDIAM(i) > (3.*PRED(i));
        CD(i)=PRED(i);
    else
        CD(i)=EQDIAM(i);
    end
end
DBH=1:size(EVT);
AGB=1:size(EVT);
AGC=1:size(EVT);
%Assumes that HT and CD are in units of meters, DBH is in units centimeters
for i=1:size(EVT);
    if EVT(i) <=3 %NaNs, Spruce-Fir Forest and Woodland (Dry Mesic), Spruce-Fir Forest and Woodland (Moist Mesic), & Forest Meadow
        DBH(i)=0.82.*(HT(i).*sqrt(CD(i)));
        AGB(i)=(HT(i).*100).*3.1415952.*((DBH(i)./2).^2).*SG(i)./1000; %Calculate aboveground biomass (AGC) in units of Kg
        AGC(i)=AGB(i).*C(i);
    elseif EVT(i) == 4 %Mixed Conifer Forest and Woodland (Dry Mesic)
        DBH(i)=0.82.*HT(i).*sqrt(CD(i));
        AGB(i)=(HT(i).*100).*pi().*((DBH(i)./2).^2).*SG(i)./1000; 
        AGC(i)=AGB(i).*C(i);
    elseif EVT(i) == 5 %Mixed Conifer Forest and Woodland (Moist Mesic)	
        DBH(i)=0.82.*HT(i).*sqrt(CD(i));
        AGB(i)=(HT(i).*100).*pi().*((DBH(i)./2).^2).*SG(i)./1000; 
        AGC(i)=AGB(i).*C(i);
    elseif EVT(i) == 7 %Blue Spruce Fringe Forest
        DBH(i)=0.82.*HT(i).*sqrt(CD(i));
        AGB(i)=(HT(i).*100).*pi().*((DBH(i)./2).^2).*SG(i)./1000; 
        AGC(i)=AGB(i).*C(i);
    elseif EVT(i) == 10 %Aspen Forest and Woodland (Dry Mesic)
        DBH(i)=0.82.*HT(i).*sqrt(CD(i));
        AGB(i)=(HT(i).*100).*pi().*((DBH(i)./2).^2).*SG(i)./1000; 
        AGC(i)=AGB(i).*C(i);
    elseif EVT(i) == 11 %Aspen Forest and Woodland (Moist Mesic)
        DBH(i)=0.82.*HT(i).*sqrt(CD(i));
        AGB(i)=(HT(i).*100).*pi().*((DBH(i)./2).^2).*SG(i)./1000; 
        AGC(i)=AGB(i).*C(i);
    elseif EVT(i) == 13 %Ponderosa Pine Forest
        DBH(i)=0.82.*HT(i).*sqrt(CD(i));
        AGB(i)=(HT(i).*100).*pi().*((DBH(i)./2).^2).*SG(i)./1000; 
        AGC(i)=AGB(i).*C(i);
    elseif EVT(i) == 14 %Gambel Oak-Mixed Montane Shrubland
        DBH(i)=0.82.*HT(i).*sqrt(CD(i));
        AGB(i)=(HT(i).*100).*pi().*((DBH(i)./2).^2).*SG(i)./1000; 
        AGC(i)=AGB(i).*C(i);
    else %Non Forested - missclassified
        DBH(i)=0.82.*HT(i).*sqrt(CD(i));
        AGB(i)=(HT(i).*100).*pi().*((DBH(i)./2).^2).*SG(i)./1000; 
        AGC(i)=AGB(i).*C(i);
    end
end
DBH=transpose(DBH);
AGC=transpose(AGC);
AGB=transpose(AGB);
CD=transpose(CD);

At this point I'm ready to export the files back into a CSV or XLSX format

I suggest using the dlmwrite command in Matlab to write the export_utm.m outputs of the file to a *.CSV

redondo=[ID,UTME,UTMN,EVT,HT,PRED,EQDIAM,MAJAX,MINAX,AREA,MAXHT,MINHT,MEANHT,CD,DBH,SG,C,AGB,AGC];
dlmwrite('J:\Work\PostDoc2\Brooks_Harpold\Jemez\QGIS\STEMS\redondo.csv',redondo,'delimiter',',','precision','%.3f');

Note: dlmwrite does not allow you to copy the text header directly.

Alternately, use the xlswrite command in Matlab to write to an XLSX file with header:

col_header={'ID','UTME','UTMN','EVT','HT','PRED','EQDIAM','MAJAX','MINAX','AREA','MAXHT','MINHT','MEANHT','CD','DBH','SG','C','AGB','AGC'};
redondo=[ID,UTME,UTMN,EVT,HT,PRED,EQDIAM,MAJAX,MINAX,AREA,MAXHT,MINHT,MEANHT,CD,DBH,SG,C,AGB,AGC];
xlswrite('D:\Work\PostDoc2\Brooks_Harpold\Jemez\QGIS\STEMS\redondo.xlsx',redondo,1,'A2');
xlswrite('D:\Work\PostDoc2\Brooks_Harpold\Jemez\QGIS\STEMS\redondo.xlsx',col_header,1,'A1');

Note: You can only do this for ~1 million rows in an XLSX file.

Calculating the total biomass on a per area basis, derived from the individual trees:

After calculating the individual tree biomass I was also interested in calculating the total biomass on a per pixel basis. I kept this analysis at a larger pixel scale, e.g. 10m and 30m, because individual trees tend to take up a significant amount of space with their crowns and their root balls.

In QGIS I used the Vector > Research Tools > Vector Grid option to generate a polygon grid that aligns with the raster topographic surface models:

After I generated the gridded polygon I used the Vector > Analysis Tools > Points in Polygon module to sample the sum of the above ground carbon (calculated from the equations above) for each given EVT/species.

The final step is to set the viewing to a graduated color ramp:

Load the Catchment polygons (Upper Jaramillo, History Grove, La Jara, and ZOB) into QGIS and use the Select by Location feature.

Create new shape files for each set of stems.

 

 

References

Bright, B. C., Hicke, J. A., & Hudak, A. T. (2012). Estimating aboveground carbon stocks of a forest affected by mountain pine beetle in Idaho using lidar and multispectral imagery. Remote Sensing of Environment, 124, 270-281. 

Erdody, T. L., & Moskal, L. M. (2010). Fusion of LiDAR and imagery for estimating forest canopy fuels. Remote Sensing of Environment, 114(4), 725-737. 

Falkowski, M. J., Smith, A. M., Gessler, P. E., Hudak, A. T., Vierling, L. A., & Evans, J. S. (2008). The influence of conifer forest canopy cover on the accuracy of two individual tree measurement algorithms using lidar data.Canadian Journal of Remote Sensing, 34(sup2), S338-S350.

Hudak, A. T., Lefsky, M. A., Cohen, W. B., & Berterretche, M. (2002). Integration of lidar and Landsat ETM+ data for estimating and mapping forest canopy height. Remote sensing of Environment, 82(2), 397-416. 

Hudak, A. T., Crookston, N. L., Evans, J. S., Hall, D. E., & Falkowski, M. J. (2008). Nearest neighbor imputation of species-level, plot-scale forest structure attributes from LiDAR data. Remote Sensing of Environment, 112(5), 2232-2245.

 Kashian, D. M., Turner, M. G., & Romme, W. H. (2005). Variability in leaf area and stemwood increment along a 300-year lodgepole pine chronosequence.Ecosystems, 8(1), 48-61. 

Lefsky, M. A., Cohen, W. B., Hudak, A., Acker, S. A., & Ohmann, J. (1999). Integration of lidar, Landsat ETM+ and forest inventory data for regional forest mapping. International Archives of Photogrammetry and Remote Sensing, 32, 119-126. 

Litton, C. M., Ryan, M. G., Tinker, D. B., & Knight, D. H. (2003). Belowground and aboveground biomass in young postfire lodgepole pine forests of contrasting tree density. Canadian Journal of Forest Research, 33(2), 351-363.

Litton, C. M., Ryan, M. G., & Knight, D. H. (2004). Effects of tree density and stand age on carbon allocation patterns in postfire lodgepole pine. Ecological Applications, 14(2), 460-475. 

Swetnam, T. L. (2013). Cordilleran forest scaling dynamics and disturbance regimes quantified by aerial LiDAR. THE UNIVERSITY OF ARIZONA. 

Swetnam, T. L., & Falk, D. A. (2014). Application of Metabolic Scaling Theory to reduce error in local maxima tree segmentation from aerial LiDAR. Forest Ecology and Management, 323, 158-167.

Turner, M. G., Tinker, D. B., Romme, W. H., Kashian, D. M., & Litton, C. M. (2004). Landscape patterns of sapling density, leaf area, and aboveground net primary production in postfire lodgepole pine forests, Yellowstone National Park (USA). Ecosystems, 7(7), 751-775.