# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Processamento de dados LiDAR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# Autor: Humberto Menecheli
#        Departamento de Ciências Florestais
#        ESALQ/USP - 09/Jun/2022
#
# Processa nuvens LiDAR (plantio/parcela) coletadas com drone em Goiás 
#
# Linguagem de programação:
#       R (v 4.2)
#       pacote lidR* (v 4.0.2 May 4th, 2022) Jean Romain Roussel
#
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
rm(list=ls(all=TRUE))                                   # Limpa memória
gc()

# Carrega pacotes lidR e terra
require(lidR); require (terra)

# Define diretório de trabalho ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lidarDir<- "C:/LiDAR/"

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ESTUDO DA NUVEM de pontos LiDAR obtidos abaixo da copa em parcela
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lidarPar <- "eparcela.las"

# Define nome completo do aquivo com a nuvem LiDAR da parcela ~~~~~~~~~
parcelArq <- paste0(lidarDir, lidarPar)

# Leitura parcial dos dados LiDAR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
par1 <- readLAS(parcelArq,              # lê apenas XYZ, e classificação
                select = "xyzc")  
print(par1@data)
summary(par1@data)

# A nuvem de pontos pertence a uma parcela circular com área de 400m2.
# Para conferir o perímetro dessa parcela, vamos plotar alguns pontos
# nessa borda. Primeiramente, podemos definir uma função que cria n
# coordenadas na borda de um círculo com centro e raio conhecidos.
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
borda <- function(x_centro, y_centro, raio, n_pontos){
  angulo <- seq(0, 2 * pi, length = n_pontos)
  x      <- raio * cos(angulo) + x_centro
  y      <- raio * sin(angulo) + y_centro
  pontos <- cbind(x, y)
  return(pontos)
}

# Cálculo do centro da parcela ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
xcentral <- round(mean(par1@data$X), digits =1)
ycentral <- round(mean(par1@data$Y), digits =1)

# Uso da função borda para gerar as coordenadas de 20 pontos no
# perímetro da parcela círcular de área igual a 400 m2,
# ou seja, raio =~ 11.18 m2
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
raio  <- 11.18
coord <- borda(xcentral, ycentral, raio, 20)

# Inclusão das coordenadas do ponto central ~~~~~~~~~~~~~~~~~~~~~~~~~~~
coord <- rbind(coord, c(xcentral, ycentral))

# A função vect(), do pacote "terra", permite criar um objeto espacial
# para as coordenadas geradas. Definiremos como CRS dessas coordenadas,
# o mesmo atribuído originalmente à nuvem de pontos (EPSG:31982).
crsNuvem <- par1@crs$input
pts <- vect(coord, type="points", crs = crsNuvem)

# As coordenadas podem agora ser plotadas em um mapa 2D do pacote terra
terra::plot(pts, col="red", pch=20, cex=2)

# A nuvem de pontos original pode ser "desbastada" (p.ex. 100 pts/m2)
# com a função decimate_points() e exibida como um plot 3D
parDesb   <- decimate_points(par1, random(100))
parDesb3D <- plot(parDesb,
                  color  = "Z",
                  bg     = "black")

# As 21 coordenadas definidas para demarcar o perímetro e o centro
# da parcela podem ser exibidas a uma certa altura no mesmo plot 3D.
# Para isso, definimos essa altura (p.ex.: 35m) como um atributo em um 
# data-frame de uma única coluna "H" com 21 linhas ...
alturas <- data.frame(H=rep(c(35), each = 21))

# ... incluímos essas alturas como atributo das coordenadas,
# e recriamos o objeto espacial que as representa no formato sf
pts <- sf::st_as_sf(vect(coord, 
                         atts=alturas, 
                         type="points", 
                         crs = crsNuvem)
                    )
# ... para, então, acrescentá-las no plot 3D
add_treetops3d(parDesb3D, pts, z="H")


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ESTUDO DA NUVEM DE PONTOS gerada por sobrevoo LiDAR acima DO PLANTIO
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lidarPla <- "eplantio.las"

# Define nome completo do aquivo com a nuvem LiDAR do plantio ~~~~~~~~~
plntioArq <- paste0(lidarDir, lidarPla)

# Leitura parcial dos dados LiDAR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pla <- readLAS(plntioArq, select = "xyzc")
print(pla@data)
summary(pla@data)

# Clipa parcela circular de 400 m2 de dentro da nuvem do plantio ~~~~~~
par2  <- clip_circle(pla, xcentral, ycentral, raio)

# Plot da parcela colorida por classificação ~~~~~~~~~~~~~~~~~~~~~~~~~~
parPlantio <- plot(par2,
                   size   = 3,
                   color  = "Classification",
                   bg     = "black")
# Inclusão das esferas delimitadoras na parcela extraída do plantio ~~~
add_treetops3d(parPlantio, pts, z="H") # H define a altura das esferas

# Plot da nuvem desbastada de pontos do plantio ~~~~~~~~~~~~~~~~~~~~~~~
plaDesb   <- decimate_points(pla, random(100))
plaDesb3D <- plot(plaDesb,
                  color  = "Z",
                  bg     = "black")
# Inclusão das esferas delimitadoras no plantio ~~~~~~~~~~~~~~~~~~~~~~~
add_treetops3d(plaDesb3D, pts, z="H") # H define a altura das esferas

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Este script foi criado em Junho/2022 depois de instalar as mais
# recentes versões dos pacotes lidR, sf, raster e terra nessa data.
#
# * lidR latest development version: sequência de instalação de pacotes
#   devtools::install_github("Jean-Romain/rlas", dependencies=TRUE)
#   install.packages("raster", dependencies = TRUE)
#   install.packages('terra', repos='https://rspatial.r-universe.dev')
#   devtools::install_github("Jean-Romain/lidR")
#
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

---

Referências sobre o pacote lidR para processamento de dados LiDAR:

• The lidR package (eBook)
• Roussel et al.(2020) lidR: An R package for analysis of ALS data
• RDocumentation (lidR)
• Vignettes

e sobre os pacotes terra e sf para processamento de dados espaciais:

• The terra package
• Simple Features for R