OVERVIEW
The purpose of this script is to calculate the Hegyi’s competition Index from individual tree outputs that come from TASS (growth model). Here, we first present the methodology broken down, and then present a loop to extract this information from all years for mean values.
TASS stands for "Tree and Stand Simulator", and it is a tree model that is biologically based and spatially explicit. There are currently three forms of TASS: (1) TASS III, which is begins to extend TASS into more complex stand structures and multi-age and multi-species cohorts; (2) TASS II, which is the current main version of the software that "predicts the potential growth and yield of even-aged, single species managed stands for ten commercial tree species"; and (3) TIPSY (Table Interpolation for Stand Yields), which provides access to yield tables generated in TASS II.
For a full list of the ten species for which TASS estimates for and for more information, visit the government of BC's webpage here: https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-inventory/growth-and-yield-modelling/tree-and-stand-simulator-tass
In this script, we will be estimating metrics for Douglas-fir, Pseudotsuga menziesii (code = Fdc).
The Hegyi’s competition Index was first established in 1974 and designed to help in the management of Jack-pine stands. The index is expressed as:
\(C_{i} = \Sigma\frac{D_{j}/D_{i}}{R_{ij}}\)
where \(i\) is the target tree and \(j\) is the competitor tree, \(D\) is the DBH (Diameter at Breast Height) of the tree \(i\) or \(j\), and \(R_{ij}\) is the distance between trees \(i\) and \(j\). Note, the original equation presented by Hegyi in 1974 added one to the denominator:
\(C_{i} = \Sigma\frac{D_{j}/D_{i}}{R_{ij} + 1}\)
(but note that in Garci 2022, it is suggested this is an additional FOOT, not METER).
Higher values suggest greater competition, whereas lower values indicate lesser competition.
The library siplab can be used to calculate this index, following protocols outlined by Garci (2022).
Hegyi, F. A simulation model for managing Jack-pine stands. In Growth Models for Tree and Stand Simulation; Fries, J., Ed.; Research Note No. 30; Royal College of Forestry, Department of Forest Yield Research: Stockholm, Sweden, 1974; pp. 74–9
Garci, O. 2022. Computing Hegyi’s (and other) Competition Indices --- siplab Vignette #1
BEGIN DATA PROCESSING
Note: if any of the libraries are not yet installed, install below using the sample code provided.
install.packages("siplab")## Warning: package 'siplab' is in use and will not be installedThe data file we will be working with is titled "interview.tre". Locate where this file is saved on your computer and modify the location in the code below as necessary prior to running the code.
Note: the TASS outputs have many lines of descriptor code prior to the table. In order to get the data in table form, we can use the "skip" function to ignore the descriptors at the top. This descriptor can have useful metadata, however, including: date data extracted; plot description; units for volume, height, diameter, basal area and crown area; top height info; predominant height info; and plot size.
For this file, the metadata is:
***** TASS v2.07.76 Custom Tree Summary ***** Tue Apr 12 10:59:19 2022 Plot Description: Units: volume: cubic m, height: m, diameter: cm, basal area: square m, crown area: square m. Top height: average height of 0 trees/ha of largest diameter. Predominant height: average height of 0 tallest trees/ha. Plot size: 49.41m x 49.41m = 0.2441ha
TASSdata<-read.table(file = "C:\\Users\\coras\\Documents\\Stored Online\\3. Work Contracts\\BC FLNRORD Stats 2\\Interview_Assignment\\Interview.tre" , header = FALSE, sep = "", dec = ".", skip = 12) # file specifies where the file is saved on the local computer; header indicates whether to include the column names or not; sep indicates how each column is separated; dec indicates what character is used to indicate decimal places; and skip is used to indicate how many lines of code to skip before reading the table data
head(TASSdata) #Note that with how we brought this data in, there are no column names## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
## 1 1 Fdc 12 4.51 4.364 1.642 0.004 0 0.873 1.123
## 2 2 Fdc 12 5.90 4.781 2.217 0.007 0 0.624 2.870
## 3 3 Fdc 12 4.75 4.317 1.825 0.004 0 0.873 4.367
## 4 4 Fdc 12 4.88 4.411 2.050 0.004 0 0.624 6.613
## 5 5 Fdc 12 4.08 3.837 1.666 0.003 0 1.123 8.360
## 6 6 Fdc 12 5.59 4.861 2.107 0.006 0 0.873 10.356Because the column names were separated on multiple rows, I opted not to bring the column names in. Instead, we will name the columns accordingly.
This is the column headings plus the first row of data from the 'Interview.tre' file:
IDSPECTOTA CROWN
CODE AGE DBH HEIGHT WIDTHVOLUMEVOLUME ROW COLUMN ------------------------------------------------------------ 1 Fdc 12 4.51 4.364 1.642 0.004 0.000 0.873 1.123
As we can see, some of the code names are jumbled together. This is the order they should be in (acquired from the output.in file): DEF CUSTOM_TREE trees HEADERS ID_CODE xxxxx SPECIES xxxx TOTAL_AGE xxxx DBH xxxx.xx HEIGHT xxx.xxx CROWN_WIDTH xxx.xxx VOLUME xx.xxx VOLUME xx.xxx TOP_DIB 10 STUMP 0.30 ROW xxx.xxx COLUMN xxx.xxx END
TASSdata<-TASSdata %>% rename(ID_CODE = V1,
SPECIES = V2,
TOTAL_AGE = V3,
DBH = V4,
HEIGHT = V5,
CROWN_WIDTH = V6,
VOLUME = V7,
VOLUME2 = V8, # TOP_DIB 10 STUMP 0.30; Note: VOLUME was already used as a column header, so changed it here
ROW = V9,
COLUMN =V10)
#Check column names updated correctly
head(TASSdata)## ID_CODE SPECIES TOTAL_AGE DBH HEIGHT CROWN_WIDTH VOLUME VOLUME2 ROW COLUMN
## 1 1 Fdc 12 4.51 4.364 1.642 0.004 0 0.873 1.123
## 2 2 Fdc 12 5.90 4.781 2.217 0.007 0 0.624 2.870
## 3 3 Fdc 12 4.75 4.317 1.825 0.004 0 0.873 4.367
## 4 4 Fdc 12 4.88 4.411 2.050 0.004 0 0.624 6.613
## 5 5 Fdc 12 4.08 3.837 1.666 0.003 0 1.123 8.360
## 6 6 Fdc 12 5.59 4.861 2.107 0.006 0 0.873 10.356str(TASSdata)## 'data.frame': 4755 obs. of 10 variables:
## $ ID_CODE : int 1 2 3 4 5 6 7 8 9 10 ...
## $ SPECIES : chr "Fdc" "Fdc" "Fdc" "Fdc" ...
## $ TOTAL_AGE : int 12 12 12 12 12 12 12 12 12 12 ...
## $ DBH : num 4.51 5.9 4.75 4.88 4.08 5.59 4.6 4.18 4.66 4.33 ...
## $ HEIGHT : num 4.36 4.78 4.32 4.41 3.84 ...
## $ CROWN_WIDTH: num 1.64 2.22 1.82 2.05 1.67 ...
## $ VOLUME : num 0.004 0.007 0.004 0.004 0.003 0.006 0.005 0.003 0.004 0.004 ...
## $ VOLUME2 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ROW : num 0.873 0.624 0.873 0.624 1.123 ...
## $ COLUMN : num 1.12 2.87 4.37 6.61 8.36 ...Inspect the data. We know that we have information at 10-year intervals, for example.
#table(TASSdata$ID_CODE) #ID code is duplicated many times because at multiple time steps
table(TASSdata$TOTAL_AGE) #We can see the distribution of number of trees as the forest gets older##
## 12 22 32 42 52 62 72 82 92 102
## 709 673 661 607 497 423 354 307 270 254We will need to separate the data by which year is being simulated in order to calculate the Hegyi’s competition Index for each time step. Separating them out here allows for easy determination of the index by any given year; however, this could also be done using a loop to pull out by year, which will be demonstrated later.
TASSdata_12<-subset(TASSdata, TASSdata$TOTAL_AGE==12)
TASSdata_22<-subset(TASSdata, TASSdata$TOTAL_AGE==22)
TASSdata_32<-subset(TASSdata, TASSdata$TOTAL_AGE==32)
TASSdata_42<-subset(TASSdata, TASSdata$TOTAL_AGE==42)
TASSdata_52<-subset(TASSdata, TASSdata$TOTAL_AGE==52)
TASSdata_62<-subset(TASSdata, TASSdata$TOTAL_AGE==62)
TASSdata_72<-subset(TASSdata, TASSdata$TOTAL_AGE==72)
TASSdata_82<-subset(TASSdata, TASSdata$TOTAL_AGE==82)
TASSdata_92<-subset(TASSdata, TASSdata$TOTAL_AGE==92)
TASSdata_102<-subset(TASSdata, TASSdata$TOTAL_AGE==102)To calculate Hegyi’s competition Index, we will follow the protocol/methodology as outlined in Garci (2022). Because the competition index is only relevant within a given timeframe, we will calculate it for each time-step separately. We will run the code for age 12 initially, then repeat for other years of interest.
We know that the information we need is the dbh of the target and competitor tree (DBH), and the distance between trees (calulated using the ROW and COLUMN information).
We first must convert the data to a "spatstat" point pattern object of class ppp. We can accomplish this using the function ppp(), and include information on the coordinate vectors (here, ROW AND COLUMN), the plot dimensions (which we obtain from the metdata, 49.41m x 49.41m), and the 'marks' (DBH).
trees_12<-ppp(TASSdata_12$ROW, TASSdata_12$COLUMN, c(0,49.41), c(0,49.41), marks = TASSdata_12$DBH)
summary(trees_12) #We see we have 709 points, with an average intensity of 0.29 points per square unit. Median value of 5, mean pf 5.04 and max of 7.9. These median, mean and max values pertain to the DBH within the data.## Marked planar point pattern: 709 points
## Average intensity 0.2904133 points per square unit
##
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
##
## marks are numeric, of type 'double'
## Summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.020 4.510 5.000 5.043 5.520 7.930
##
## Window: rectangle = [0, 49.41] x [0, 49.41] units
## Window area = 2441.35 square unitsNow that we have information in individual trees, we want to determine information on the sizes and distances for pairs of trees (e.g., "pairwise"). In siplab, we will use the function "kernel". One consideration is what radius around each tree for which competition is relevant needs to be specified. sometimes, maxR = Rmax, or other times it is specified according to a desired number of nearest neighbours (maxN = n). Within siplab, we can use the powers_ker() function, but we must choose and specify the R (radius) to be considered. Garci (2022) use 6 m. Here, we will provide an example using 6 m.
hegyi_12 <- pairwise(trees_12, maxR=6, kernel=powers_ker,
kerpar = list(pi=1, pj=1, pr=1, smark=1)) #These parameters happen to be the default, so the kerpar function is not necessarily needed here. It is included in case changing parameters is desired.
head(marks(hegyi_12))## mark cindex
## 1 4.51 3.437300
## 2 5.90 3.352014
## 3 4.75 5.402397
## 4 4.88 5.427710
## 5 4.08 7.486211
## 6 5.59 5.159041summary(hegyi_12 ) #709 points, mean cindex = 8.495## Marked planar point pattern: 709 points
## Average intensity 0.2904133 points per square unit
##
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
##
## Mark variables: mark, cindex
## Summary:
## mark cindex
## Min. :3.020 Min. : 2.650
## 1st Qu.:4.510 1st Qu.: 7.276
## Median :5.000 Median : 8.489
## Mean :5.043 Mean : 8.495
## 3rd Qu.:5.520 3rd Qu.: 9.793
## Max. :7.930 Max. :17.178
##
## Window: rectangle = [0, 49.41] x [0, 49.41] units
## Window area = 2441.35 square unitsIf it is desired to use the originally published equation that has the additional 1 foot (note again, here we use meters), then the following code can be used:
#hegyi_12b <- pairwise(trees_12, maxR=6, kernel=hegyiorig_ker) #No kerpar needed; Note: as of April 18, 2022, "hegyiorig_ker" is not found...
#head(marks(hegyi_12b))Another consideration is that of edge effects. Because trees on the edges of the plots may experience lesser competition in the index because they have fewer neighbours, they may wish to be excluded from the calculation. Because we used 6 m as our max R, we will also use this value as the minimum distance from the edge for a tree to be included in the measurement.
hegyi_12_trim<-edges(hegyi_12, -6)
summary(hegyi_12_trim) # down from 709 to 428 pointa. Mean cindex = 9.396## Marked planar point pattern: 428 points
## Average intensity 0.3058217 points per square unit
##
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
##
## Mark variables: mark, cindex
## Summary:
## mark cindex
## Min. :3.020 Min. : 5.578
## 1st Qu.:4.487 1st Qu.: 8.278
## Median :4.945 Median : 9.236
## Mean :5.026 Mean : 9.396
## 3rd Qu.:5.485 3rd Qu.:10.421
## Max. :7.800 Max. :17.178
##
## Window: rectangle = [6, 43.41] x [6, 43.41] units
## (37.41 x 37.41 units)
## Window area = 1399.51 square unitsThe above procedure can be repeated for any year within the data. For example, let's compare the c-index from year 12 (above) to year 102.
#First get data into point pattern object
trees_102<-ppp(TASSdata_102$ROW, TASSdata_102$COLUMN, c(0,49.41), c(0,49.41), marks = TASSdata_102$DBH)
#Next, get the pairwise values, with cindex calculated, using the same radius for comparison
hegyi_102 <- pairwise(trees_102, maxR=6, kernel=powers_ker,
kerpar = list(pi=1, pj=1, pr=1, smark=1))
head(marks(hegyi_102))## mark cindex
## 1 17.64 3.1574838
## 2 42.95 1.0635969
## 3 34.29 1.5144719
## 4 30.50 1.4116411
## 5 21.51 1.8145853
## 6 42.98 0.7461124summary(hegyi_102 ) #254 points, mean cindex = 3.3339## Marked planar point pattern: 254 points
## Average intensity 0.1040409 points per square unit
##
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
##
## Mark variables: mark, cindex
## Summary:
## mark cindex
## Min. : 9.25 Min. :0.6213
## 1st Qu.:17.03 1st Qu.:1.7222
## Median :23.20 Median :2.9199
## Mean :26.40 Mean :3.3339
## 3rd Qu.:35.59 3rd Qu.:4.6156
## Max. :62.79 Max. :9.3315
##
## Window: rectangle = [0, 49.41] x [0, 49.41] units
## Window area = 2441.35 square units#Remove the trees within 6 m of the edge
hegyi_102_trim<-edges(hegyi_102, -6)
summary(hegyi_102_trim) #152 points, mean cindex 3.5754## Marked planar point pattern: 152 points
## Average intensity 0.1086096 points per square unit
##
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
##
## Mark variables: mark, cindex
## Summary:
## mark cindex
## Min. : 9.25 Min. :0.6729
## 1st Qu.:17.09 1st Qu.:1.9597
## Median :23.45 Median :3.3325
## Mean :26.58 Mean :3.5754
## 3rd Qu.:33.85 3rd Qu.:5.0221
## Max. :62.79 Max. :9.0291
##
## Window: rectangle = [6, 43.41] x [6, 43.41] units
## (37.41 x 37.41 units)
## Window area = 1399.51 square unitsIf you wish to create a table for getting the mean and median c-index each year, you can create a loop and insert into a table. This first table does not remove the trees closest to the edge.
unique(TASSdata$TOTAL_AGE)## [1] 12 22 32 42 52 62 72 82 92 102AGES<-c("12", "22", "32", "42", "52", "62", "72", "82", "92", "102")
#Create table for outputs to be put into
table.hegyi <- data.frame (matrix (ncol = 6, nrow = 0))
colnames (table.hegyi) <- c ("AGES", "N", "Mean_CIndex", "Median_CIndex", "Min_CIndex", "Max_CIndex")
for (i in 1:length(AGES)) {
temp.data<-TASSdata %>% dplyr::filter(TOTAL_AGE == AGES[i])
temp.output_b<-ppp(temp.data$ROW, temp.data$COLUMN, c(0,49.41), c(0,49.41), marks = temp.data$DBH)
temp.output <- pairwise(temp.output_b, maxR=6, kernel=powers_ker,
kerpar = list(pi=1, pj=1, pr=1, smark=1))
#Extract values desired
x0<- AGES[i]
x1<- temp.output$n
x2<- mean(temp.output$marks$cindex)
x3<- median(temp.output$marks$cindex)
x4<- min(temp.output$marks$cindex)
x5<- max(temp.output$marks$cindex)
tab.hegyi.temp <- cbind.data.frame(AGES=x0, N=x1, Mean_CIndex =x2, Median_CIndex = x3, Min_CIndex = x4, Max_CIndex = x5)
table.hegyi<-rbind(table.hegyi, tab.hegyi.temp)
}Inspect output.
table.hegyi## AGES N Mean_CIndex Median_CIndex Min_CIndex Max_CIndex
## 1 12 709 8.494628 8.489255 2.6501458 17.17800
## 2 22 673 8.244051 7.965267 2.6539459 27.17890
## 3 32 661 8.552913 7.918938 2.3790064 21.48527
## 4 42 607 8.099468 7.632372 1.9470197 22.91535
## 5 52 497 6.636577 6.314706 1.3706437 21.67530
## 6 62 423 5.670498 5.350773 1.0674957 16.00414
## 7 72 354 4.743097 4.376623 0.8814882 13.90842
## 8 82 307 4.109454 3.702253 0.7498891 12.60243
## 9 92 270 3.554745 3.125575 0.6365349 10.54228
## 10 102 254 3.333916 2.919871 0.6212693 9.33146Repeat the above while accounting for edge effects.
unique(TASSdata$TOTAL_AGE)## [1] 12 22 32 42 52 62 72 82 92 102AGES<-c("12", "22", "32", "42", "52", "62", "72", "82", "92", "102")
#Create table for outputs to be put into
table.hegyi.trim <- data.frame (matrix (ncol = 6, nrow = 0))
colnames (table.hegyi.trim) <- c ("AGES", "N", "Mean_CIndex", "Median_CIndex", "Min_CIndex", "Max_CIndex")
for (i in 1:length(AGES)) {
temp.data<-TASSdata %>% dplyr::filter(TOTAL_AGE == AGES[i])
temp.output_b<-ppp(temp.data$ROW, temp.data$COLUMN, c(0,49.41), c(0,49.41), marks = temp.data$DBH)
temp.output_c <- pairwise(temp.output_b, maxR=6, kernel=powers_ker,
kerpar = list(pi=1, pj=1, pr=1, smark=1))
temp.output<-edges(temp.output_c, -6)
#Extract values desired
x0<- AGES[i]
x1<- temp.output$n
x2<- mean(temp.output$marks$cindex)
x3<- median(temp.output$marks$cindex)
x4<- min(temp.output$marks$cindex)
x5<- max(temp.output$marks$cindex)
tab.hegyi.temp <- cbind.data.frame(AGES=x0, N=x1, Mean_CIndex =x2, Median_CIndex = x3, Min_CIndex = x4, Max_CIndex = x5)
table.hegyi.trim<-rbind(table.hegyi.trim, tab.hegyi.temp)
}Inspect output.
table.hegyi.trim## AGES N Mean_CIndex Median_CIndex Min_CIndex Max_CIndex
## 1 12 428 9.396471 9.235927 5.5780290 17.178000
## 2 22 398 9.013589 8.535864 4.6272104 27.178895
## 3 32 389 9.321679 8.498133 3.3215933 21.485269
## 4 42 360 8.874755 8.484074 2.4298357 22.915346
## 5 52 294 7.213252 6.960952 1.7831489 21.675297
## 6 62 251 6.157237 6.074541 1.3184394 16.004137
## 7 72 218 5.315127 5.202661 1.0473446 13.908416
## 8 82 186 4.517093 4.216686 0.8325629 12.602431
## 9 92 162 3.836386 3.495725 0.7023294 10.542277
## 10 102 152 3.575380 3.332486 0.6729121 9.029086Here, we calulated the Hegyi Competition Index from a TASS output for Douglas Fir within a 43.41 x 43.41 m plot at 10-year intervals from year 12 to 102. We first broke the process down for how to inspect for each year separately, and then created a loop to inspect all years at once for easier comparison. Some important considerations that need to be made are (1) what radius is being specified for being included in the calculation between pairs of trees and (2) whether edge effects are to be accounted for, removing the trees within the determined radius from (1) from the forest edge. Both results broadly demonstrate that the Hegyi Competition Index reduced as forest age increases, as does the number of living trees.