Lab08 – Coarse Woody Debris (volume of residual wood after harvest)

Lab08 Data Download

Lab 08: Coarse woody debris volume remaining after harvest

[Excerpt from Brian Davis’ 2017 thesis,”REFINEMENT OF A DRONE-BASED METHOD FOR ESTIMATING COARSE WOODY DEBRIS AND BIOMASS RESIDUE FOLLOWING FOREST HARVEST”]

"… Pine forests of the southern US are a major source of woody debris residues used for bioenergy (Milbrandt 2005). “In 2012, approximately 40% of the world’s biomass fuel, in the form of wood pellets created from CWD harvest, was exported from the United State to Europe (REN21 2013), most from the southern US. Estimates are that 60% of the world renewable energy could be in the form of biomass fuel by 2030 (IRENA 2014). With such a large and growing market for the biomass fuel, the need for accurate, quick estimation of CWD residue has never been greater. From an operational standpoint, contractors need means of estimating biomass to plan harvests and schedule equipment. From a regulatory standpoint, estimates of retained biomass are needed to ensure compliance with biomass retention guidelines…"

How to use data from drone flight to determine the volume of residual wood left after harvest???

For this lab,I want to know the volume (cubic-meters) of residual wood within the Lab8_ClpBndy shapefile.  Upload a Word document with the volume.

 

library("lidR")
library("rgdal")
library("EBImage")
# shapepath<- "U:/FANR5640/Spring2019/labs/Lab08Data"
# shapename<- "Lab8_ClpBndy"
# shapef<- file.path(shapepath,shapename)
# shape_spdf<- readOGR(dsn = shapepath, layer = shapename)
# shape_p<- shape_spdf@polygons[[1]]@Polygons[[1]]
# 
# 
# mypath<- "U:/FANR5640/Spring2018/FANR5640/lab9a"
# myfile<- "exhibitb.las"
# infile<- file.path(mypath,myfile)
# lastmp = readLAS(infile, select = "xyzirc")##, Intensity = F, ReturnNumber = F, NumberOfReturn = F, ScanAngle = F, EdgeOfFlightline = F, ScanDirectionFlag = F, UserData = F, PointSourceId = F, color = T)
# las = lasclip(lastmp, shape_p)
# las_original = las
# plot(las)
# summary(las)
# plot(las)
# writeLAS(las,file.path(shapepath,"clplas.las"))

#####Load shapefile boundary
shapepath<- "U:/FANR5640/Spring2019/labs/Lab08Data"
shapename<- "PileBoundary"
shapef<- file.path(shapepath,shapename)
shape_spdf<- readOGR(dsn = shapepath, layer = shapename)
shape_p<- shape_spdf@polygons[[1]]@Polygons[[1]]

#####Read and clip LAS to shapefile boundary
myfile<- "clplas.las"
infile<- file.path(shapepath,myfile)
lastmp = readLAS(infile, select = "xyzirc")##, Intensity = F, ReturnNumber = F, NumberOfReturn = F, ScanAngle = F, EdgeOfFlightline = F, ScanDirectionFlag = F, UserData = F, PointSourceId = F, color = T)
las = lasclip(lastmp, shape_p)
las_original = las
plot(las)
summary(las)

#####Classify Ground
plot(las, color = "Classification")
las<- las_original
las<- lasground(las,csf(cloth_resolution = 0.5, class_threshold = 0.25), last_returns = FALSE)
las_ground <- las
plot(las, color = "Classification")


#####NORMALIZE
las_original<- las
las <- lasnormalize(las, knnidw(k = 8, p = 2))
las_normal<- las
summary(las_original)
summary(las_normal)
plot(las_original)
plot(las_normal, color = "Classification")

#####Remove below-ground points
las_normal_nolow<- lasfilter(las_normal,(Z >= 0))
plot(las_normal_nolow)
summary(las_normal_nolow)

#####Create canopy height model; run following commands with res = 2.5, 1, 0.5, 0.1
chm <- grid_canopy(las_normal_nolow, res = 0.1, pitfree(c(0,2,5,10,15), c(0, 1.5)))
plot(chm)
writeRaster(chm,filename = file.path(shapepath,"chm_01m.tif"), format="GTiff",overwrite=TRUE)

At this point, you have saved the canopy height model out as a raster layer you can open in ArcMap.  Use the Spatial Analyst’s Raster Calculator to (possibly) refine or remove low-lying noise.  For your last calculation, use the Raster Calculator to convert the heights into cubic volumes.  The pixel values of your raster CHM represent the height above the ground in meters.  If the layer has a 1-meter cell size, then the volume of each pixel would be calculated by

width * Length * height or (1 * 1 * <pixel value>)

Look at the layer properties to determine the pixel size.  Once you have made your calculation, then, use the ZONAL STATISTICS AS TABLE command to summarize the volume for the area within the PileBoundary layer.