Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- lidar:projetos:rcode_invdrone [2022/06/13 13:27] – criada lcer
+++ lidar:projetos:rcode_invdrone [2024/03/23 20:17] (atual) – edição externa 127.0.0.1
@@ Linha 1: / Linha 1: @@
+[[ lidar:projetos |{{ :lidar:projetos:retornar.png?nolink&25x25 }}]]
+<code>
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# 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")
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+</code>
+----
+Referências sobre o pacote //lidR// para processamento de dados LiDAR:
+  * [[https://r-lidar.github.io/lidRbook/ | The lidR package (eBook)]]
+  * [[https://hal.inrae.fr/hal-02972284/file/roussel_2020.pdf | Roussel et al.(2020) lidR: An R package for analysis of ALS data]]
+  * [[https://www.rdocumentation.org/packages/lidR/versions/4.0.1 | RDocumentation (lidR)]]
+  * [[https://cran.rstudio.com/web/packages/lidR/vignettes/ | Vignettes]]
+e sobre os pacotes //terra// e //sf// para processamento de dados espaciais:
+  * [[https://rspatial.org/terra/pkg/index.html | The terra package]]
+  * [[https://r-spatial.github.io/sf/articles/sf1.html | Simple Features for R ]]