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)
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')