11 Sensoriamento remoto

11.1 Imagens de RPAs

Este capítulo oferece uma introdução ao processamento e classificação de imagens de RPAs no ambiente R, usando algoritmos de aprendizado capazes de realizer diversos processamentos. Ele também fornece um tutorial de referência conciso e prático, que oferece aos leitores uma nosso geral do que é possível realizar no Sistema R com o processamento e classificação de imagens de RPAs.

É necessário instalar alguns pacotes necessários para aplicar esta tarefa:

install.packages("devtools", dependency=TRUE)
library(devtools)
install_github("filipematias23/FIELDimageR")
install.packages("sp", dependency=TRUE)
install.packages("raster", dependency=TRUE)
install.packages("rgdal", dependency=TRUE)

Carregar os pacote que foram baixados:

library(FIELDimageR)
library(raster)

Baixar a imagem ortomosaicada:

EX1<-stack("https://www.dropbox.com/s/qxg6fn9ljy35ssm/EX1_RGB.tif?dl=1")

Plotar a imagem nas bandas RGB:

plotRGB(EX1, r = 1, g = 2, b = 3)

Remover o solo e trabalhar apenas com a vegetação para aplicar os índices de vegetação RGB:

EX1.RemSoil<- fieldMask(mosaic = EX1, Red = 1, Green = 2, Blue = 3, index = "HUE")

Aplicado os índices de vegetação. Aplicaremos o índice NGRDI, BGI e podemos criar um índice usando as bandas disponíveis. Criaremos como exemplo myIndex com a fórmula Red-Blue/Green:

EX1.Indices<- indices(mosaic = EX1.RemSoil$newMosaic, Red = 1, Green = 2, Blue = 3, 
                      index = c("NGRDI","BGI"), myIndex = c("(Red-Blue)/Green"))

11.2 Imagens de satélites

Esse procedimento será realizado com imagens de satélite (Sentinel 2), porêm pode ser aplicado com imagens de RPA, desde que sejam multiespectrais.

Carregar pacotes necessários para trabalhar com os dados raster. Caso não tenha algum dos pacotes, realize a sua instalação:

library(raster)
library(knitr)
library(sp)
library(rgdal)
library(ggplot2)
library(viridis)
library(rasterVis)
library(LSRS)

Baixar o arquivo . sentinel2.tif.

imagesentinel=raster('D:/livro/TudodoRa/sentinel2.tif')

imagesentinel=raster('sentinel2.tif')

Visualizar os dados:

imagesentinel

É necessário criar camadas individuais para cada uma das bandas espectrais:

b1 <- raster('sentinel2.tif', band=1)
b2 <- raster('sentinel2.tif', band=2)
b3 <- raster('sentinel2.tif', band=3)
b4 <- raster('sentinel2.tif', band=4)
b5 <- raster('sentinel2.tif', band=5)
b6 <- raster('sentinel2.tif', band=6)
b7 <- raster('sentinel2.tif', band=7)
b8 <- raster('sentinel2.tif', band=8)
b9 <- raster('sentinel2.tif', band=9)
b10 <- raster('sentinel2.tif', band=10)
b11 <- raster('sentinel2.tif', band=11)
b12 <- raster('sentinel2.tif', band=12)

Comparar duas bandas para ver se elas possuem a mesma extensão:

compareRaster(b2, b3)

Plotar a banda 4 para pré-visualização:

plot(b4)

image(b4)

Visualizar a imagem nas bandas do RGB:

RGB <- stack(list(b4, b3, b2))              
plotRGB(RGB, axes = TRUE, stretch = "lin", main = "Sentinel RGB colour composite")

Juntar todas as bandas num só arquivo:

t <- stack(b1,b2, b3, b4, b5, b6, b7, b8, b9, b10, b11, b12)
st <- brick('sentinel2.tif')
plot(st)

Aplicar o índice de vegetação NDVI, para o Sentinel 2 com: NIR = 8, red = 4.

Criar a VI (vegetation index) por meio de função:

VI <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}

NDVI:

ndvi <- VI(st, 8, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Sentinel2-NDVI")

Outras fórmula de aplicar o NDVI

vi2 <- function(x, y) {
  (x - y) / (x + y)
}
ndvi2 <- overlay(st[[8]], st[[4]], fun=vi2)
plot(ndvi2, col=rev(terrain.colors(10)), main="Sentinel2-NDVI")

Visualizar o NDVI em histograma

hist(ndvi,
     main = "Distribuição dos valores de NDVI",
     xlab = "NDVI",
     ylab= "Frequência",
     col = "wheat",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

Visualizar apenas a vegetação com NDVI acima de 0.4:

veg <- reclassify(ndvi, cbind(-Inf, 0.4, NA))
plot(veg, main='Vegetação')

Reclassificar o NDVI e difini-lo por classes numéricas:

vegc <- reclassify(ndvi, c(-Inf,0.25,1, 0.25,0.3,2, 0.3,0.4,3, 0.4,0.5,4, 0.5,Inf, 5))
plot(vegc,col = rev(terrain.colors(4)), main = 'NDVI reclassificado')

Criar uma classificação não supervisionada a partir do NDVI:

Converter o raster (NDVI) a um vetor/matriz:

nr <-getValues(ndvi)
str(nr)

É importante definir o gerador de pontos, porque o “kmeans” inicia os centros em locais aleatórios:

set.seed(99)

Criar 10 clusters, permitir 500 iterações, comece com 5 conjuntos aleatórios usando o método Lloyd:

kmncluster <- kmeans(na.omit(nr), centers = 10, iter.max = 500,
                     nstart = 5, algorithm = "Lloyd")

Ver o vetor/matriz:

str(kmncluster)

Crie uma cópia do NDVI para não perder os dados:

knr <- ndvi

Agora substitua os valores das células de varredura pelo kmncluster$cluster:

knr[] <- kmncluster$cluster

Realize o plot do NDVI e do kmeans:

par(mfrow = c(1, 2))
plot(ndvi, col = rev(terrain.colors(10)), main = "NDVI")
plot(knr, main = "Kmeans", col = viridis_pal(option = "D")(10))

Se quiser traçar a classificação kmeans ao lado da renderização do RGB para verificar a qualidade da classificação e identificação das classes:

par(mfrow = c(1, 2))
plotRGB(RGB, axes = FALSE, stretch = "lin", main = "RGB")
plot(knr, main = "Kmeans", yaxt = 'n', col = viridis_pal(option = "D")(10))

Aplicar outros índices de vegetação com o pacote LSRS:

NDVI=NDVI(b8,b4)
SAVI=SAVI(b8,b4)
TGSI=TGSI(b4,b2,b3)
MSAVI=MSAVI(b8,b4, Pixel.Depth=1)
EVI=EVI(b8,b4,b2,Pixel.Depth=1)
NBR=NBR(b8,b11)

par(mfrow = c(3, 2))
plot(NDVI,lwd=4,main="NDVI",xlab="easting", ylab="northing")
plot(SAVI,lwd=4,main="SAVI",xlab="easting", ylab="northing")
plot(TGSI,lwd=4,main="TGSI",xlab="easting", ylab="northing")
plot(MSAVI,lwd=4,main="MSAVI",xlab="easting", ylab="northing")
plot(EVI,lwd=4,main="EVI",xlab="easting", ylab="northing")
plot(NBR,lwd=4,main="NBR",xlab="easting", ylab="northing")

11.3 Curvas de nível e modelo 3D a partir do Modelo Digital de Elevação

Primeiro carregar os pacotes necessários:

library(raster)
library(plot3D)

Carregar o dado raster do pacote Raster (Volcano) para ser usado como exemplo:

filled.contour(volcano, color.palette = terrain.colors)

Criar as curvas de nível:

cont <- contourLines(volcano)
fun <- function(x) x$level
LEVS <- sort(unique(unlist(lapply(cont, fun))))
COLS <- terrain.colors(length(LEVS))

Plotar somente as curvas de nível:

contour(volcano)

Plotar o modelo 3D com curvas de nível:

x <- seq(1, nrow(volcano), by = 3)
y <- seq(1, ncol(volcano), by = 3) 
Volcano <- volcano [x, y]

Exemplo 1:

ribbon3D(z = Volcano, contour = TRUE, zlim= c(-100, 200),image = TRUE)

Exemplo 2:

persp3D(z = Volcano, contour = TRUE, zlim= c(-200, 200), image = FALSE)

Exemplo 3:

persp3D(z = Volcano, x = x, y = y, scale = FALSE,contour = list(nlevels = 20, col = "red"),zlim = c(-200, 200), expand = 0.2,image = list(col = grey (seq(0, 1, length.out = 100))))

Exemplo 4:

persp3D(z = Volcano, contour = list(side = c("zmin", "z", "350")), zlim = c(-100, 400), phi = 20, image = list(side = 350))

Exemplo 5:

persp3D(z = volcano, shade = 0.3, col = gg.col(100))

11.4 LIDAR

O exemplo a seguir é um processamento de imagens LIDAR, por meio da qual será segmentado árvores individuais e métricas.

Carregar pacote:

library(lidR)
library(raster)

Baixar dados Example.las

Carregar uma área florestal de exemplo a partir de imagens LIDAR que será trabalhada:

las = readLAS('C:/Users/Jefferson/Dropbox/Livro/Example.las')
plot(las)

O lasground fornece vários algoritmos para classificar os pontos de referência. Essa função é conveniente para gráficos de pequeno a médio porte, como o que estamos processando.

Classificar pontos do solo (pontos de referência):

las = lasground(las, csf())
plot(las, color = "Classification")

É necessário definir o terreno em 0 metros. Deve-se subtrair o MDT para obter pontos de aterramento em 0, mas aqui não será usado um MDT, mas vamos interpolar exatamente cada ponto.

Definir altura normalizada:

las = lasnormalize(las, tin())
plot(las)

Na próxima etapa, será usado um algoritmo que requer um modelo de altura da copa a partir da nuvem de pontos.

Calcular um modelo de altura das copas:

algo = pitfree(thresholds = c(0,10,20,30,40,50), subcircle = 0.2)
chm  = grid_canopy(las, 0.5, algo)

Plotar o CHM:

plot(chm, col = height.colors(50))

A segmentação pode ser alcançada com lastrees. Aqui foi escolhido o algorítmo de bacia com um limiar de 4 metros. A nuvem de pontos foi atualizada e cada ponto agora tem um número que se refere a uma árvore individual (treeID). Pontos que não são árvores recebem o valor de ID NA.

Realizar a segmentação das árvores:

algo = watershed(chm, th = 4)
las  = lastrees(las, algo)

Remova os pontos que não estão atribuídos a uma árvore:

trees = lasfilter(las, !is.na(treeID))

Plotar a segmentação:

plot(trees, color = "treeID", colorPalette = pastel.colors(100))

Calcular algumas métricas:

hulls  = tree_hulls(las, func = .stdmetrics)
spplot(hulls, "zmax")

No exemplo anterior, mesmo se a segmentação for feita usando um modelo de altura do dossel, a classificação foi feita na nuvem de pontos. Isso ocorre porque lidR é uma biblioteca orientada à nuvem de pontos. Mas pode-se querer que o raster trabalhe com rasters. Nesse caso, a função divisor de águas pode ser usada de forma independente:

crowns = watershed(chm, th = 4)()
plot(crowns, col = pastel.colors(100))

Criar poligonos de contornos a partir da copa:

contour = rasterToPolygons(crowns, dissolve = TRUE)

Plotar o CHM e contornos:

plot(chm, col = height.colors(50))
plot(contour, add = T)

11.5 Modelo Digital de Elevação MDE

Carregar pacotes:

library(raster)

Carregar dados e plotar:

MDE=raster("https://www.dropbox.com/s/b2rzimq500rmj5o/MDE.tif?dl=1")
plot(MDE)

Usar a função terrain() para calcular/extrair algumas informações topográficas:

Declividade <- terrain(MDE, "slope")
Aspecto <- terrain(MDE, "aspect")
TPI <- terrain(MDE, "TPI") # Topographic Position Index (Índice de posição topográfica)
TRI <- terrain(MDE, "TRI") # Terrain Ruggedness Index (Índice de robustez do terreno)
Rugosidade <- terrain(MDE, "roughness")
Escoamento <- terrain(MDE, "flowdir")
Hillshade <- hillShade(Declividade, Aspecto, angle=45, direction=0, filename='', normalize=FALSE)

Juntar todos os dados com a função stack():

topo <- stack(MDE, Declividade, Aspecto, TPI, TRI, Rugosidade, Escoamento, Hillshade)

Renomear os dados para aparecerem no plot:

names(topo) <- c("MDE", "Declividade","Aspecto", "TPI", "TRI", "Rugosidade", "Escoamento", "Hillshade")

Plotar os dados:

plot(topo)