In tankwin08/waveformlidar: Waveform LiDAR Data Processing and Analysis

stopifnot(require(knitr))
options(width = 90)
opts_chunk$set(
  comment = NA,
  message = FALSE,
  warning = FALSE,
  eval = if (isTRUE(exists("params"))) params$EVAL else FALSE,
  dev = "png",
  dpi = 150,
  fig.asp = 0.8,
  fig.width = 5,
  out.width = "60%",
  fig.align = "center"
)
library(waveformlidar)
#ggplot2::theme_set(theme_default())

Goal

The purpose of this vignette is to show an example to conduct waveform decompostion using differnt models such as the Gaussian, adaptive Gaussian and Weibull functions algorithms using waveformlidar package.

A Case Study

As a case study, we will use the waveform LiDAR data of the waveformlidar package, which comes from the National Ecological Observatory Networks (NEON). The example dataset used here is a subset (500 waveforms) from the Harvard Forest, Massachusetts, USA. Each row represents one waveform and each column represents a time bin. The temporal resolution of the time bin can be 1/2/4 nanosecond (ns), which is up to the data collecting system. In this case, the temporal resolution is 1 ns with the vertical distance approximately 0.15m.

data(return, package = "waveformlidar")
head(return)

Overview

The deconvolution is an algorithm-based process that is to reverse the effect of convolution on the recorded signals, and the decomposition is a process which can provide estimates of the location and properties of objects along the pulse (Wagner et al. 2006). In this package, two deconvolution methods including the Gold and Richardson-Lucy (RL) are available for conducting the deconvolution. The detailed description of these algorithms can be found in (Zhou et al. 2017b). For the decomposition, both the outgoing pulse and the return pulse are nearly following some probability distribution such as Gaussian distribution in terms of shape. Thus, the information inherent in waveform can be extracted through fitting waveforms with a mixture of models with the specific distribution or waveform components. By interpreting these models' parameters or waveform components, the targets such as vegetation and ground interacting with outgoing pulse along the path can be characterized.

Prepare data and decompostion

library(data.table)
data(return)
wf<-data.table(index=c(1:nrow(return)),return)
###decomposition, we chose two examples: one is the simple waveform and another was more complex waveform
r1<- decom(wf[1,])    ##default, with smooth was applied
r2<- decom(wf[1,], smooth = FALSE)  ##use the raw waveform

##for more complicated waveform
r3<- decom(wf[182,])
###when the waveform is mixed with too much noise, we should consider use more filtering options
r4<- decom(wf[182,],smooth=TRUE,width=5)

###we also can fit the waveform with the other models such as adaptive Gaussian and Weibull
r5<-decom.adaptive(wf[182,])

r6<-decom.weibull(wf[182,])

Decomposition results

To visually check the decomposition results using differnt models, we plot two sample waveform(one is simple waveform, one is comlicated waveform). The plot was generated using the plot form the base. The ggplot2 is also another option for demonstrating results.

###prepare the data for three functions
#########fig 3a
y0<- return[182,]
y0[y0==0]<-NA
y1<- as.numeric(y0 - min(y0,na.rm = T) +1)
y1<-y1[1:wavelen(y1)]

re<-decom(wf[182,],smooth=TRUE, width=5)[[2]]
yi<-cbind(re[1:3,2],re[4:6,2],re[7:9,2])

sumyi<-0
x<-1:wavelen(y1)
for (i in 1:3){
  sumyi<-sumyi + yi[i,1] * exp(-(x - yi[i,2])**2/(2 * yi[i,3]**2))
}


####fig 3b.
ayi<-r5[[3]]

sumayi<-0
x<-1:wavelen(y1)
for (i in 1:nrow(ayi)){
  sumayi<-sumayi + ayi[i,2] * exp(-abs(x - ayi[i,3])**ayi[i,5]/(2 * ayi[i,4]**2))
}


#############fig. 3c
wyi<-r6[[3]]

sumwyi<-0
x<-1:wavelen(y1)
for (i in 1:nrow(wyi)){
  sumwyi<-sumwyi + wyi[i,2] * wyi[i,5]/wyi[i,4] * abs(x - wyi[i,3])^(wyi[i,5] - 1)/(wyi[i,4])^(wyi[i,5] - 1) * exp(-abs(x - wyi[i,3])^wyi[i,5]/(wyi[i,4])^wyi[i,5])
}


###############begin to plot
setwd("A:/research/2018/waveform_summary_R/waveformlidar/vignette") ##change to where you want to save your results
par(family = 'serif')
png('differnt_decomposition_results.png', width = 9.25*2.5, height = 5.73*1.2, units = 'in',res = 500)

par(mfrow=c(1,3), oma=c(0,0,0,0),mai=c(1.0,1.0,0.8,0))
###fig3.a

x<-1:wavelen(y1)
plot(seq_along(y1),y1,type="l",col="red",xlab="(a) Time(ns)",ylab="Intensity",
     cex.lab=2.5,lwd=3.5,axes=F,main="Gaussian",cex.main=2,ylim=c(-30,230))
axis(side=1, at=seq(0,(wavelen(y1)+15),15),cex.axis=2.5)
axis(side=2, at=seq(-30,(max(y1)+40),20),cex.axis=2.5)
lines(seq_along(y1),y1,type="l",col="red",lwd=3)

###idea waveform components
for (i in 1:3){
  lines(x,yi[i,1] * exp(-(x - yi[i,2])**2/(2 * yi[i,3]**2)),lty=2,col=i+2,lwd=3)
}

lines(x,sumyi,lty=3,col="black",lwd=3)
legend("topright",legend=c("1","2","3"),col=c(3,4,5),lty=c(2,2,2),lwd=3,cex=2,box.lwd="white")
legend("topleft",legend=c("SW","RW"),col=c("black","red"),lty=c(2,1),lwd=3,cex=2,box.lwd="white")

#####fig. 3b
plot(seq_along(y1),y1,type="l",col="red",xlab="(b) Time(ns)",ylab="Intensity",
     cex.lab=2.5,lwd=3.5,axes=F,main="Adaptive Gaussian",cex.main=2,ylim=c(-30,max(sumayi)))
axis(side=1, at=seq(0,(wavelen(y1)+15),15),cex.axis=2.5)
axis(side=2, at=seq(0,(max(sumayi)+20),20),cex.axis=2.5)
lines(seq_along(y1),y1,type="l",col="red",lwd=3)

###idea waveform components
for (i in 1:nrow(ayi)){
  lines(x,ayi[i,2] * exp(-abs(x - ayi[i,3])**ayi[i,5]/(2 * ayi[i,4]**2)),lty=2,col=i+2,lwd=3)
}

lines(x,sumayi,lty=3,col="black",lwd=3)
legend("topright",legend=c("1","2","3"),col=c(3,4,5),lty=c(2,2,2),lwd=3,cex=2,box.lwd="white")
legend("topleft",legend=c("SW","RW"),col=c("black","red"),lty=c(2,1),lwd=3,cex=2,box.lwd="white")


#############################fig.3c weibull function


plot(seq_along(y1),y1,type="l",col="red",xlab="(c) Time(ns)",ylab="Intensity",
     cex.lab=2.5,lwd=3.5,axes=F,main="Weibull",cex.main=2,ylim=c(-30,230))
axis(side=1, at=seq(0,(wavelen(y1)+15),15),cex.axis=2.5)
axis(side=2, at=seq(0,(max(y1)+40),20),cex.axis=2.5)
lines(seq_along(y1),y1,type="l",col="red",lwd=3)


for (i in 1:nrow(wyi)){
  lines(x,wyi[i,2] * wyi[i,5]/wyi[i,4] * (x - wyi[i,3])^(wyi[i,5] - 1)/(wyi[i,4])^(wyi[i,5] - 1) * exp(-(x - wyi[i,3])^wyi[i,5]/(wyi[i,4])^wyi[i,5]),
        lty=2,col=i + 2,lwd=3)
  #lines(x,intens[i] * exp(-abs(x - ti[i])^2/(2 * std[i]^2)),lty=2,col="gray",lwd=0.5)
}
lines(x,sumwyi,lty=3,col="black",lwd=3)
legend("topright",legend=c("1","2","3"),col=c(3,4,5),lty=c(2,2,2),lwd=3,cex=2,box.lwd="white")
legend("topleft",legend=c("SW","RW"),col=c("black","red"),lty=c(2,1),lwd=3,cex=2,box.lwd="white")


dev.off()

Waveform processing over an region

The above example just show the waveform processing over an individual waveform, while most of time we need to process millions of waveforms for a region. The following show you an example how to implent it with apply function.

dr3<-apply(wf,1,decom.adaptive)

####to collect all data
rfit3<-do.call("rbind",lapply(dr3,"[[",1)) ## waveform is correctly decomposed with index,some are not correct index by NA
ga3<-do.call("rbind",lapply(dr3,"[[",2))   ###the original results, which can help to get more detailed results.
pa3<-do.call("rbind",lapply(dr3,"[[",3))   ###useful decompostion results for next step or geolocation transformation.

####delete some wrong ones
rid<-rfit3[!is.na(rfit3),]
wid<-setdiff(as.numeric(unlist(wf[,1])),rid)  ###index of waveforms needs to be reprocessed
rpars<-pa3[!is.na(pa3[,1]),]