lasR: Simple lidar reading, processing, viewing

Documented in experimental_metrics

#'@title
#'  compute experimental lidar metrics
#'
#'@description
#'  take advantage of xy and 3d information
#'
#'@details
#'  <Delete and Replace>
#'
#'\cr
#'Revision History
#' \tabular{ll}{
#'1.0 \tab 4/4/2019 \cr
#'}
#'
#'@author
#'
#'Jacob Strunk <Jstrunk@@fs.fed.us>
#'
#'@param x lidR object X vector : las1$X
#'@param y lidR object Y vector : las1$Y
#'@param z lidR object Z vector : las1$Z
#'@param i lidR object intensity vector : las1$Intensity
#'@param rn lidR object return vector : las1$ReturnNumber
#'@param cl lidR object return vector : las1$Classification
#'@param r lidR object red vector : las1$R
#'@param g lidR object green vector : las1$G
#'@param b lidR object blue vector : las1$B
#'@param resxy resolution to use when computing xy statistic - grid resolution
#'@param resxyz resolution to use when computing xyz statistic - voxel resolution
#'@param ressurf resolution to use when computing raster areas - grid resolution
#'@param htcover height to use when computing cover and threshold for eliminating ground returns
#'@param zthresh height thresholds to use for height ratios round(100*length(z[ z > zthresh[i]]) / length(z[ z > htcover]),2)
#'@param adjustz (optional) method used to shift z values up/down. When there is no DTM or the DTM is poor (e.g. when used with DAP) this can help
#'@param outlier a two value vector with minimum and maximum values used to restrict the range of z values considered
#'
#'@return
#'  list of statistics, currently only computed on x,y,z
#'
#'@examples
#'    if(!"las1" %in% ls()) las1 = readLAS("D:\\Box\\sync\\R\\analyses\\wa_dsm_env_fia\\data\\plot_clips_ht\\plot_44.laz")
#'    test2 = lasmetrics(las1,experimental_metrics(X,Y,Z))
#'    test2
#'
#@import some_package,some_package2
#'
#'@export
#
#'@seealso \code{\link{another_function}}\cr \code{\link{yet_another_function}}\cr

#Desired upgrades to this function:
#
#

# x = function(x){}

#copy function arguments and use this code to format arguments
##writeClipboard(paste(gsub("^[[:space:]]*[,]*","#'@param ",gsub("=.*"," ",readClipboard())),collapse="\n"))


experimental_metrics = function(
  x=NA
  ,y=NA
  ,z=NA
  ,i=NA
  ,rn=NA
  ,cl=NA
  ,r=NA
  ,g=NA
  ,b=NA
  ,resxy=20
  ,resxyz = 20
  ,ressurf=5
  ,htcover = 6
  ,zthresh = c(3,6,12,20)
  ,adjustz = c("none","positive","min0","zq01")
  ,outlier=c(min=NA,max=NA)
  ){

  #adjust z values to e.g. set the minimum value to zero
  #especially useful with DAP over a coarse ground model...
  if(!is.na(z[1]))if(!is.null(adjustz[1]))if(!is.na(adjustz[1])){
    if(adjustz[1]=="positive") if(min(z) < 0) z = z - min(z) + 1
    if(adjustz[1]=="min0") z = z - min(z)+.001
    if(adjustz[1]=="zq01") z = z - quantile(z,.01)
  }

  if(!is.na(outlier[1])){
    outlier_ids = z < outlier[1] | z > outlier[2]
    z = z[!outlier_ids]
    if(!is.na(x[1])) x = x[!outlier_ids]
    if(!is.na(y[1])) y = y[!outlier_ids]
    if(!is.na(i[1])) i = i[!outlier_ids]
    if(!is.na(rn[1])) rn = rn[!outlier_ids]
    if(!is.na(cl[1])) cl = cl[!outlier_ids]
    if(!is.na(r[1])) r = r[!outlier_ids]
    if(!is.na(g[1])) g = g[!outlier_ids]
    if(!is.na(b[1])) b = b[!outlier_ids]
  }

  mets_in = list()

  if(!is.na(x[1]) & !is.na(y[1]) & !is.na(z[1])){

    library(plyr)
    library(raster)

    x1 = x[z >= htcover]
    y1 = y[z >= htcover]
    z1 = z[z >= htcover]

    rtxy = (x*y)^(1/2)
    rtxy1 = (x1*y1)^(1/2)
    rtxyz = (x*y*z)^(1/3)
    rtxyz1 = (x1*y1*z1)^(1/3)
    xyz=data.frame(x = round(x/resxy) *  resxy, y = round(y/resxy) * resxy, z = z)
    xyz1=data.frame(x = round(x1/resxy) *  resxy, y = round(y1/resxy) * resxy, z = z1)
    combnsxy = count(data.frame(round(x/resxy), round(y/resxy)))
    combnsxy1 = count(data.frame(round(x1/resxy), round(y1/resxyz)))
    combnsxyz = count(data.frame(round(y/resxyz), round(x/resxyz), round(x/resxyz)))
    combnsxyz1 = count(data.frame(round(y1/resxyz), round(x1/resxyz), round(z1/resxyz)))

    r1=raster(res=ressurf,xmn=floor(min(x)),xmx=floor(min(x)) + ceiling((max(x)-min(x))/ressurf)*ressurf,ymn=floor(min(y)),ymx=floor(min(y)) + ceiling((max(y)-min(y))/ressurf)*ressurf)
    ch=rasterize(cbind(x,y),r1,field=z,fun=function(x,...)quantile(x,.975,na.rm=T))
    chspdf = as(ch,"SpatialPixelsDataFrame")

    #typical z metrics
    mets_in["zMin"] = min(z1,na.rm=T)
    mets_in["zMinRat"] = round(100*mets_in[["zMin"]] / min(z,na.rm=T) ,2)
    mets_in["zMax"] = max(z1,na.rm=T)
    mets_in["zMaxRat"] = round(100*mets_in[["zMax"]] / max(z,na.rm=T) , 2)
    mets_in["zSd"] = sd(z1,na.rm=T)
    mets_in["zSdRat"] = round(100*mets_in[["zSd"]] / sd(z,na.rm=T) ,2)
    mets_in["zCv"] = cv(z1,na.rm=T)
    mets_in["zCvRat"] = round(100*mets_in[["zCv"]] / cv(z,na.rm=T) ,2)
    mets_in["zIQ"] = diff(range(z1,na.rm=T))
    mets_in["zIQRat"] = round(100*mets_in[["zIQ"]]/ diff(range(z,na.rm=T)),2)
    mets_in["zCover"] = round(100*length(z1)/ length(z),2)
    mets_in["zCoverHt"] = htcover
    mets_in[paste("zq",c(1,5,10,20,30,40,50,60,70,80,90,95,99),sep="")] = quantile(z1,c(1,5,10,20,30,40,50,60,70,80,90,95,99)/100)
    mets_in[paste("zRat",zthresh,sep="")] = lapply(zthresh,function(thr,z1)round(100*sum(z>thr)/length(z1),2),z1)

    mets_in["xMin"] = min(x1,na.rm=T)
    mets_in["xMax"] = max(x1,na.rm=T)
    mets_in["xSd"] = sd(x1,na.rm=T)
    mets_in["xRatHt"] = round(100*diff(range(x1,na.rm=T)) / diff(range(x,na.rm=T)),2)

    mets_in["yMin"] = min(y1,na.rm=T)
    mets_in["yMax"] = max(y1,na.rm=T)
    mets_in["ySd"] = sd(y1,na.rm=T)
    mets_in["yRatHt"] = round(100*diff(range(y1,na.rm=T)) / diff(range(y,na.rm=T)),2)

    mets_in["xyMinRt"] = min(rtxy1,na.rm=T)
    mets_in["xyMinRtRat"] = round(100*mets_in[["xyMinRt"]] / min(rtxy,na.rm=T),2)
    mets_in["xyMaxRt"] = max(rtxy1,na.rm=T)
    mets_in["xyMaxRtRat"] = round(100*mets_in[["xyMaxRt"]] / max(rtxy,na.rm=T),2)
    mets_in["xySdRt"] = sd(rtxy1,na.rm=T)
    mets_in["xySdRtRat"] = round(100*mets_in[["xySdRt"]] / sd(rtxy,na.rm=T),2)
    mets_in["xyCvRt"] = cv(rtxy1,na.rm=T)
    mets_in["xyCvRtRat"] = round(100*mets_in[["xyCvRt"]] / cv(rtxy,na.rm=T),2)
    mets_in["xyCor"] = cor(x1,y1)
    mets_in["xyCorRat"] = round(100*mets_in[["xyCor"]] / cor(x,y),2)
    mets_in["xyArea"] = nrow(combnsxy1) * resxy^2

    mets_in["xyzMinRt"] = min(rtxyz1)
    mets_in["xyzMinRtRat"] = round(100*mets_in[["xyzMinRt"]] / min(rtxyz,na.rm=T),2)
    mets_in["xyzMaxRt"] = max(rtxyz1)
    mets_in["xyzMaxRtRat"] = round(100*mets_in[["xyzMaxRt"]] / max(rtxyz,na.rm=T),2)
    mets_in["xyzSdRt"] = sd(rtxyz1)
    mets_in["xyzSdRtRat"] = round(100*mets_in[["xyzSdRt"]] / sd(rtxyz,na.rm=T),2)
    mets_in["xyzCvRt"] = cv(rtxyz1)
    mets_in["xyzCvRtRat"] = round(100*mets_in[["xyzCvRt"]] / cv(rtxyz,na.rm=T),2)
    mets_in["xyzCor"] = cor(rtxy1,z1)
    mets_in["xyzCorRat"] = round(100*mets_in[["xyzCor"]] / cor(rtxy,z),2)

    mets_in["voxVol"] = nrow(combnsxyz1)* resxyz^3
    mets_in["voxVolRat"] = round(100*mets_in[["voxVol"]] / (nrow(combnsxyz)* resxyz^3),2)
    #compute volume between min / max for given raster cell:
    mets_in["gridVol"] = sum(aggregate(z~x+y,data=xyz1,FUN = function(z,resxy)(max(z,na.rm=T) - min(z,na.rm=T))*resxy^2,resxy=resxy)[,3])
    mets_in["gridVolRat"] = round(100*mets_in[["gridVol"]] / sum(aggregate(z~x+y,data=xyz,FUN = function(z,resxy)(max(z,na.rm=T) - min(z,na.rm=T))*resxy^2,resxy=resxy)[,3]),2)
    mets_in["areaCover"] = round(100*mets_in[["xyArea"]] / ( nrow(combnsxy) * resxy^2 ),2)
    mets_in["surfArea"] = surfaceArea(chspdf)
    mets_in["surfAreaRat"] = round(mets_in[["surfArea"]] / (sum(ch[]>0,na.rm=T)*ressurf^2)*100,2)

  }

  #other ways to compute metrics
  if(F){
    mets_in["xyarea"] = length(unique( round(y/resxy) + round(x/resxy) / 5000 )) * resxy^2
    require(raster)
    r1=raster(resxy=resxy,xmn=floor(min(x)),xmx=floor(min(x)) + ceiling((max(x)-min(x))/resxy)*resxy,ymn=floor(min(y)),ymx=floor(min(y)) + ceiling((max(y)-min(y))/resxy)*resxy)
    mets_in["xyarea"] = length(unique(cellFromXY(r1,cbind(x,y))))*resxy*resxy

  }

  return(mets_in)

}

# all_metrics = function(obj){
#
#   cbind(
#     lasmetrics(obj,.stdmetrics)
#     ,lasmetrics(obj,entropy(Z-min(Z)+1))
#     ,lasmetrics(obj,.LAD)
#     )
#
# }



if(F) {
  require(lidR)
  if(!"las1" %in% ls()) lidR::las1 = readLAS("D:\\Box\\sync\\R\\analyses\\wa_dsm_env_fia\\data\\plot_clips_ht\\plot_66399.laz")
  test2 = lasmetrics(las1,experimental_metrics(X,Y,Z,adjustz="none",outlier=c(-15,450)))
  test2
}

if(F){

  library(lidR)
  test=list()
  test["test1"] = 1
  test
  las1 = readLAS("D:\\Box\\sync\\R\\analyses\\wa_dsm_env_fia\\data\\plot_clips_ht\\plot_44.laz")

  test1 = lasmetrics(las1,.stdmetrics)
  test2 = lasmetrics(las1,.stdmetrics)
  test2 = lasmetrics(las1,LAD(Z))
  test3 = all_metrics(las1)

  test2 = lasmetrics(las1,experimental_metrics(X,Y,Z))


  m1 = grid_metrics(las, stdmetrics(X,Y,Z,Intensity,ReturnNumber,Classification,dz=1))

}