Import Packages

library(lidR)
## Loading required package: raster
## Loading required package: sp
library(raster)
library(purrr)

Import Polygon Plot Shapefile

centShape <-shapefile("/Users/lavran_pagano/Downloads/CreeksheadPlotPolygon2.shp")
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum D_Unknown_based_on_GRS80_ellipsoid in CRS definition:
## +proj=utm +zone=17 +ellps=GRS80 +units=m +no_defs

Import list of canopy height models and mosaic them into one large raster covering the study area

#path to chm folder
mypath<- '/Users/lavran_pagano/Downloads/Lidar_Training_Data/Creeks_Head/Creeks_HeadCHMs'
#make a list of chm files
chmlist<-list.files(mypath,full.names = T, pattern = '.tif$') %>%  
map(raster)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
#mosaic raster
mosaic<- do.call(merge,chmlist)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

Make sure both the raster and the shapefile are in the same projection. In our case we will be working in Michiagn State Plane

#define projection (state plane)
statePlane<-"+proj=lcc +lat_0=41.5 +lon_0=-84.3666666666667 +lat_1=43.6666666666667
+lat_2=42.1 +x_0=3999999.999984 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=ft +no_defs"
#change chm projection
crs(mosaic) <- statePlane
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
#change shapefile projection
centShape<-spTransform(centShape,statePlane)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

Create list of pixel values

#create list of pixel values
pixelvalues<-raster::extract(mosaic,centShape,method= 'simple')
#Omit NAs for simplicity
pixelvalues <-lapply(pixelvalues,na.omit)

Create a dataframe of predictor varibles

Height covariets

#max hight
Height_Max<-lapply(pixelvalues,max)
Height_Max<-unlist(Height_Max)
#min height
Height_Min<-lapply(pixelvalues,min)
Height_Min<-unlist(Height_Min)
#height 10th
percent10<-function(x){
  quantile(x, probs = c(0.1))
}
Height_10th<-lapply(pixelvalues,percent10)
Height_10th<-unlist(Height_10th)
#height 20th
percent20<-function(x){
  quantile(x, probs = c(0.2))
}
Height_20th<-lapply(pixelvalues,percent20)
Height_20th<-unlist(Height_20th)
#height 30th
percent30<-function(x){
  quantile(x, probs = c(0.3))
}
Height_30th<-lapply(pixelvalues,percent30)
Height_30th<-unlist(Height_30th)
#height 40th
percent40<-function(x){
  quantile(x, probs = c(0.4))
}
Height_40th<-lapply(pixelvalues,percent40)
Height_40th<-unlist(Height_40th)
#height 50th
percent50<-function(x){
  quantile(x, probs = c(0.5))
}
Height_50th<-lapply(pixelvalues,percent50)
Height_50th<-unlist(Height_50th)
#height 60th
percent60<-function(x){
  quantile(x, probs = c(0.6))
}
Height_60th<-lapply(pixelvalues,percent60)
Height_60th<-unlist(Height_60th)
#height 70th
percent70<-function(x){
  quantile(x, probs = c(0.7))
}
Height_70th<-lapply(pixelvalues,percent70)
Height_70th<-unlist(Height_70th)
#height 80th
percent80<-function(x){
  quantile(x, probs = c(0.8))
}
Height_80th<-lapply(pixelvalues,percent80)
Height_80th<-unlist(Height_80th)
#height 90th
percent90<-function(x){
  quantile(x, probs = c(0.9))
}
Height_90th<-lapply(pixelvalues,percent90)
Height_90th<-unlist(Height_90th)

Crown gemoetric voulmne(cgv) cgv is calculated from the sumation of all pixel values in a plot with the minium height value subtracted from said pixel values. That number is then multiplyed by the width and length of a given pixel cell to convert to cubic meters.Our spatial resolution is 0.25ft which is equalivlent to 0.075m, Therfore given our spatail resolution cgv is calculated by sum(y-min(y)) 0.075m X 0.075m.

#CGV FUll
CGV100<-function(y){
  sum(y-min(y))* 0.005625 #0.075m X 0.075m
}
cgv100th<-lapply(pixelvalues,CGV100)
cgv100th<-unlist(cgv100th)
#CGV50th
CGV50<-function(x, y){
  y<-ifelse(x< percent50(x),x,percent50(x))
  sum(y-min(x))* 0.005625
}
cgv50th<-lapply(pixelvalues,CGV50)
cgv50th<-unlist(cgv50th)
#CGV60th
CGV60<-function(x, y){
  y<-ifelse(x< percent60(x),x,percent60(x))
  sum(y-min(x))* 0.005625
}
cgv60th<-lapply(pixelvalues,CGV60)
cgv60th<-unlist(cgv60th)
#CGV70th
CGV70<-function(x, y){
  y<-ifelse(x< percent70(x),x,percent70(x))
  sum(y-min(x))* 0.005625
}
cgv70th<-lapply(pixelvalues,CGV70)
cgv70th<-unlist(cgv70th)

Create Training Data Data Frame

#Make trainging data table
Training_data<-data.frame(Height_Max = Height_Max,
                          Height_90th = Height_90th,
                          Height_80th = Height_80th,
                          Height_70th= Height_70th,
                          Height_60th= Height_60th,
                          Height_50th=Height_50th,
                          Height_40th=Height_40th,
                          Height_30th=Height_30th,
                          Height_20th=Height_20th,
                          Height_10th=Height_10th,
                          Height_Min=Height_Min,
                          cgv100=cgv100th,
                          cgv70th= cgv70th,
                          cgv60th=cgv60th,
                          cgv50th=cgv50th)
print(Training_data)
##    Height_Max Height_90th Height_80th Height_70th Height_60th Height_50th
## 1    27.80031    23.56222    21.79986    20.18584    18.97082    17.73886
## 2    24.62012    21.42417    20.66889    19.80933    18.82830    17.53657
## 3    27.47493    23.80521    22.32861    20.66984    19.22180    17.71481
## 4    28.23649    27.48256    26.73834    26.00536    25.39520    24.42430
## 5    27.01072    23.43036    21.92004    20.49803    19.18102    18.05239
## 6    26.30650    25.28161    24.71931    24.07024    23.30112    22.34006
## 7    26.93870    23.51748    22.23682    20.89402    19.18120    17.46070
## 8    27.75986    26.73548    25.77905    24.46670    22.01528    19.53391
## 9    26.58427    23.09475    20.79948    19.04852    17.19953    15.04714
## 10   23.35016    22.23329    21.31740    20.47105    19.56552    18.49903
## 11   26.65982    24.65948    23.47535    22.45945    21.36185    20.26412
## 12   22.49252    15.62290    14.28309    12.88420    11.68808    10.61186
##    Height_40th Height_30th Height_20th Height_10th Height_Min    cgv100
## 1    16.466005    15.05236   13.781883   12.056976  6.8275189 1056.3106
## 2    15.659693    13.57484   11.807636    9.767604  4.2919598 1168.2781
## 3    15.180408    12.09247    6.960082    2.713511  0.1695230 1472.8248
## 4    23.145861    21.41590   19.601473   17.078501  9.1489496 1354.8128
## 5    17.021688    16.01139   15.022661   13.334068  6.4415989 1137.6225
## 6    21.064750    19.56804   18.127714   16.197815 12.1149921  903.2357
## 7    15.798422    14.32973   12.641495   10.280078  5.2876234 1156.2043
## 8    17.680515    16.34399   14.238638   10.138782  0.8337740 1782.2192
## 9    13.318019    11.98930   10.634428    8.677357  3.8739929 1122.9611
## 10   17.332066    16.04747   14.782765   13.229115  1.4208474 1585.5718
## 11   19.005355    17.48892   14.915391   10.766751  4.8475218 1379.1689
## 12    9.340387     7.94660    5.890117    3.224434  0.5054377  934.4627
##      cgv70th   cgv60th   cgv50th
## 1   975.7794  934.6076  880.8803
## 2  1127.1454 1094.0396 1037.3681
## 3  1401.2103 1352.4749 1285.9462
## 4  1322.9187 1302.4253 1259.7809
## 5  1067.4380 1023.1864  973.8719
## 6   876.2740  850.1863  807.8999
## 7  1095.2993 1036.6017  961.8761
## 8  1732.0441 1647.8233 1539.7400
## 9  1032.2209  969.9463  876.9298
## 10 1546.7045 1516.0229 1469.5224
## 11 1329.9808 1292.6802 1244.8591
## 12  862.7703  822.8267  775.5243

Write Training Data to a csv

write.csv(Training_data,'Creeks_Head_Covariates2.csv')