• Início
  • Conhecimento básico
  • Ferramentas
  • Criando Mapas Interativos
  • Geoprocessamento
  • Módulo IV - Geoprocessamento.

    Parte 1 - R e Geoprocessamento

    Nesta breve introdução ao geoprocessamento vamos falar de dois pacotes básicos que iremos usar com R para carregar e manipular imagens. Estes pacotes são o raster e o rgdal. Instale os dois em R usando:

    > install.packages('raster')
    > install.packages('rgdal')
    

    Vamos também carregar as seguintes imagens Sentinel 2 na resolução de 20 metros por pixel, originalmente as bandas 2, 3, 4 e 8 tem resolução de 10 metros mas foram ajustadas para 20 metros que é a resolução das demais bandas (5, 6, 7 ,8a, 11 e 12). Também baixe o modelo de elevação na resolução de 30 metros.

  • Banda 2
  • Banda 3
  • Banda 4
  • Banda 5
  • Banda 6
  • Banda 7
  • Banda 8
  • Banda 8A
  • Banda 11
  • Banda 12
  • Modelo de Elevação Digital (DEM)
  • Parte 2 - Imagens Multiespectrais

    A banda 4 corresponde ao comprimento de onda do vermelho, a 3 do verde e a 2 do azul. Ao criarmos uma imagem RGB (red, green, blue) com essas respectivas bandas, nós criamos o que chamamos de composição de cor verdadeira (true color composition - tcc em inglês).

    Caso prescise de maiores fundamentos em R, as apostilas listadas ao lado podem te ajudar.

    > library(raster)
    Loading required package: sp
    > tcc<-stack('b4lr.tif','b3lr.tif','b2lr.tif')
    > plotRGB(tcc,stretch='lin')
    

    Uma composição vnir (very near infrared) é usada para avaliar vegetação. Executamos uma composição de cor false (fcc em inglês) com as bandas 7, 6 e 5. Note que a vegetação saudável é mostrada em amarelo.

    > vnir<-stack('b7.tif','b6.tif','b5.tif')
    > plotRGB(vnir,stretch='lin')
    

    A faixa espectral correspondente ao infravermelho de ondas curtas (swir na sigla em inglês) e bastante sensível para mudanças de composição de solo e rochas. Note a separação entre as principais unidades geológicas com a composição abaixo usando as bandas 12, 11 e 8 nos canais RGB respectivamente.

    > swir<-stack('b12.tif','b11.tif','b8lr.tif')
    > plotRGB(swir,stretch='lin')
    

    Essas composições coloridas RGB são alguns exemplos. Experimente diferentes composições de bandas e compare os resultados.

    Vamos agora examinar um Modelo de Elevação Digital e os produtos que podemos extrair dele.

    Parte 3 - DEM

    Um modelo de elevação digital ou DEM é uma imagem raster onde cada pixel corresponde a um valor de elevação. O dado de estamos usando é da agência espacial japonesa JAXA do programa "ALOS World 3D - 30m (AW3D30)" com 30 metros de resolução de pixel.

    > dem<-raster('dem.tif')
    > plot(dem,col=grey(0:100/100),main='Modelo de elevação Digital - DEM')
    

    Além da imagem acima, podemos gerar a partir do DEM diversos produtos com a ajuda de R, vamos ver quais eles são aqui:

    Slope e Aspect

    A inclinação é o plano de tangência em um determinado ponto da superfície e é dividido em dois parâmetros:

  • Gradiente ou derivada ou queda: Aqui referido como 'Slope' é o ângulo que a superficie faz com o plano horizontal em radianos de 0 a pi/2.
  • Aspecto: Aqui referido como 'aspect' é a direção do megulho máximo deste plano em radianos de 0 a 2*pi em relação ao norte.
    > gradient<-terrain(dem,opt='slope',unit='radians')
    > aspect<-terrain(dem,opt='aspect',unit='radians')
    > plot(gradient,zlim=c(0,1.6),axes=FALSE)
    > plot(aspect,zlim=c(0,6.3),axes=FALSE)
    

    Slope

    Aspect

    TPI

    TPI ou índice de posição topográfica é a diferença entre o valor desta célula com a média das células a seu redor. Valores positivos indicam que essa célula é mais alta , negativo mais baixo e zero próximo de zero indicam uma superficie com um gradiente constante.

    > tpi<-terrain(dem,opt='TPI')
    > plot(tpi,zlim=c(-12,12),axes=FALSE)
    

    TRI

    Terrain Ruggedness Index (TRI) é a média da diferença absoluta entre os valores de uma célula e o valor das oito células (ou pixels) adjacentes a ele.

    > tri<-terrain(dem,opt='TRI')
    > plot(tri,zlim=c(0,30),axes=FALSE)
    

    Roughness

    Rughnees mostra a rugosidade do terreno. É expresso pela diferença dos valores máximo e do mínimo dos 8 pixels adjacentes

    > roughness<-terrain(dem,opt='roughness')
    > plot(roughness,zlim=c(0,55),axes=FALSE) 
    

    Flow

    O fluxo (flow) retorna a direção do fluxo. Ou seja, a direção de maior queda na elevação (ou o menor inclinação se todos os valores vizinhos forem mais elevados).

    > flow<-terrain(dem,opt='flowdir',unit='radians')
    > plot(flow,zlim=c(0,20),axes=FALSE)
    

  • O hillshade mostra o relevo com sombreamento, ele e gerado usando o slope e o aspecto criados acima .

    > hillshade<-hillShade(gradient,aspect,angle=45,direction=270,normalize= TRUE)
    > plot(hillshade,col=grey(0:100/100),zlim=c(0,250))
    

    Sobrepondo o hillshade na composição swir para resaltar o relevo.

    > plotRGB(swir,stretch='lin')
    > plot(hillshade,add=TRUE,legend=FALSE,axes=FALSE,col=grey(0.4:0.9),alpha=0.3,zlim=c(100,200))
    

    Você pode gravar esses dados criados acima no forma de um arquivo geotiff.

    > library(rgdal)
    > writeRaster(hillshade, filename='hillshade.tif', drive='GeoTiff')
    

    Parte 4 - PCA

    O PCA ou análise de componentes principais é um procedimento estatístico que usa uma transformação ortogonal para converter um conjuto de observações de variáveis possívelmente correlacionáveis em um conjunto de varíaveis linearmente não correlacionáveis chamado de componentes principais.

    Se existem n observações de p variáveis então o número de distintos componentes princiais será o menor valor entre n-1 e p. Esta transformação é definida de tal maneira que o primeiro componente principal tenha a maior variância possível e cada componente seguinte tenha a maior variância respeitando que seja ortogonal ao componente anterior.

    Primeiro vamos criar um raster com as bandas 2, 3, 4, 5, 6, 7, 8, 11, e 12 onde faremos a análise de componentes principais;

    > library(raster)
    Loading required package: sp
    > todas<-stack('b2lr.tif','b3lr.tif','b4lr.tif','b5.tif','b6.tif','b7.tif','b8lr.tif','b11.tif','b12.tif')
    > todas
    class      : RasterStack 
    dimensions : 5490, 5490, 30140100, 9  (nrow, ncol, ncell, nlayers)
    resolution : 20, 20  (x, y)
    extent     : 499980, 609780, 7990240, 8100040  (xmin, xmax, ymin, ymax)
    crs        : +proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
    names      :  b2lr,  b3lr,  b4lr,    b5,    b6,    b7,  b8lr,   b11,   b12 
    min values :     0,     0,     0,     0,     0,     0,     0,     0,     0 
    max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535 
    

    Vamos agora criar uma amostra do nosso conjunto de dados (população) aleatoriamente com um tamanho que seja representativo desta população. Lembre-se que análises estatísticas envolvendo imagens tem um volume de dados muito grande e consequentemente demoram um pouco (25 a 30 minutos) para serem processados, seja paciente!

    > set.seed(1)
    > amostra<-sampleRandom(todas,300000)
    

    Agora vamos efetuar o PCA na amostra para termos os preditores de cada componente e aplicar estes na nossa população gerando os nove componentes principais de nossa polulação:

    > pca <- prcomp(amostra, scale = TRUE)
    > pca
    Standard deviations (1, .., p=9):
    [1] 2.4938825 1.4405068 0.6506558 0.3245173 0.2507924 0.2101730 0.1725472
    [8] 0.1521670 0.1297261
    
    Rotation (n x k) = (9 x 9):
                PC1         PC2          PC3          PC4         PC5        PC6
    b2lr -0.3500839  0.21786006 -0.487521294  0.451152030 -0.37757283 -0.2949587
    b3lr -0.3703183  0.14580034 -0.445177787 -0.009218293  0.00742522  0.2270376
    b4lr -0.3709968  0.18873698 -0.240328170 -0.331789692  0.65117918  0.2191660
    b5   -0.3833114  0.07602207  0.178701902 -0.636503618 -0.27668479 -0.5143922
    b6   -0.3020342 -0.44426207  0.045971733 -0.193578987 -0.30085933  0.2728383
    b7   -0.2546075 -0.53012978  0.017193159  0.049618293 -0.11197494  0.2017823
    b8lr -0.2447077 -0.53589076  0.008766872  0.319041692  0.39225323 -0.3215784
    b11  -0.3505529  0.21457965  0.520360445  0.306977293  0.24457257 -0.2678492
    b12  -0.3419363  0.27972942  0.448538596  0.212598578 -0.20015435  0.5084364
                 PC7        PC8           PC9
    b2lr  0.33783895 -0.1951474  0.0835675215
    b3lr -0.69691336  0.2935101 -0.1410700518
    b4lr  0.37005967 -0.2149661  0.0606220352
    b5    0.06198892  0.2171109 -0.1344925163
    b6   -0.12450615 -0.2888241  0.6387034513
    b7    0.10786017 -0.2715668 -0.7160607841
    b8lr  0.08512160  0.5058650  0.1750840968
    b11  -0.37475020 -0.4405098 -0.0035602696
    b12   0.29081065  0.4186952 -0.0004243326
    > pci1<-predict(todas,pca,index=1)
    > pci2<-predict(todas,pca,index=2)
    > pci3<-predict(todas,pca,index=3)
    > pci4<-predict(todas,pca,index=4)
    > pci5<-predict(todas,pca,index=5)
    > pci6<-predict(todas,pca,index=6)
    > pci7<-predict(todas,pca,index=7)
    > pci8<-predict(todas,pca,index=8)
    > pci9<-predict(todas,pca,index=9)
    

    E agora vamos plotar cada um dos componentes principais criados usando:

    > par(mfrow=c(3,3))
    > plot(pci1,axes=FALSE,col=grey(0:100/100))
    > plot(pci2,axes=FALSE,col=grey(0:100/100))
    > plot(pci3,axes=FALSE,col=grey(0:100/100))
    > plot(pci4,axes=FALSE,col=grey(0:100/100))
    > plot(pci5,axes=FALSE,col=grey(0:100/100))
    > plot(pci6,axes=FALSE,col=grey(0:100/100))
    > plot(pci7,axes=FALSE,col=grey(0:100/100))
    > plot(pci8,axes=FALSE,col=grey(0:100/100))
    > plot(pci9,axes=FALSE,col=grey(0:100/100))
    

    Vamos gravar cada um destes componentes como imagem raster para uso posterior usando:

    > rm(amostra)
    > library(rgdal)
    > writeRaster(pci1,filename='pca1.tif',driver='Geotiff')
    > writeRaster(pci2,filename='pca2.tif',driver='Geotiff')
    > writeRaster(pci3,filename='pca3.tif',driver='Geotiff')
    > writeRaster(pci4,filename='pca4.tif',driver='Geotiff')
    > writeRaster(pci5,filename='pca5.tif',driver='Geotiff')
    > writeRaster(pci6,filename='pca6.tif',driver='Geotiff')
    > writeRaster(pci7,filename='pca7.tif',driver='Geotiff')
    > writeRaster(pci8,filename='pca8.tif',driver='Geotiff')
    > writeRaster(pci9,filename='pca9.tif',driver='Geotiff')
    

    Parte 5 - IHS

    Vamos agora criar imagens do tipo IHS (intensidade, matiz e saturação) a partir de uma composição RGB. Este tipo de procedimento é bastante usado para gerar uma maior resposta espectral em alguns caso, principalmente nas imagens H e S.

    Vamos gerar uma composição RGB usando SWIR (bandas 12,11 e 8) nos canais RGB inicialmente e proceder com a geração de imagens RGB a partir dela.

    > library(raster)
    > swir<-stack('b12.tif','b11.tif','b8lr.tif')
    

    Precisaremos de normalizar estas imagens para valores entre 0 e 1, fazemos isso usando:

    > Nor<-swir/65535
    > Nor
    class      : RasterBrick 
    dimensions : 5490, 5490, 30140100, 3  (nrow, ncol, ncell, nlayers)
    resolution : 20, 20  (x, y)
    extent     : 499980, 609780, 7990240, 8100040  (xmin, xmax, ymin, ymax)
    crs        : +proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
    source     : /tmp/RtmpGjvJwB/raster/r_tmp_2019-10-16_102611_32256_56357.grd 
    names      :        b12,        b11,       b8lr 
    min values :          0,          0,          0 
    max values : 0.38136873, 0.19211109, 0.08664073 
    

    Convertemos de RGB par IHS usando as seguintes fórmulas:

    I=(R+G+B)/3
    S=1-(3/(R+G+B))*mínimo(R,G,B)
    H*=acos(0.5*((R-G))+(R-B)))/sqrt((R-B)²+((R-G)*(G-B)) ) ))
    > I<-(Nor[[1]]+Nor[[2]]+Nor[[3]])/3
    > S<-1-(3/(Nor[[1]]+Nor[[2]]+Nor[[3]]))*min(Nor[[1]],Nor[[2]],Nor[[3]])
    > H<- acos((0.5*((Nor[[1]]-Nor[[2]])+(Nor[[1]]-Nor[[3]])))/(sqrt((Nor[[1]]-Nor[[2]])^2+((Nor[[1]]-Nor[[3]])*(Nor[[2]]-Nor[[3]])))))
    

    Plotando as imagens I, H e S separadamente

    > par(mfrow=c(1,3))
    > plot(I,axes=FALSE,legend=FALSE,col=grey(0:100/100))
    > plot(H,axes=FALSE,legend=FALSE,col=grey(0:100/100))
    > plot(S,axes=FALSE,legend=FALSE,col=grey(0:100/100))
    

    Criando uma Composição IHS

    > ihs<-brick(I,H,S)
    > plotRGB(ihs,stretch='lin')
    

    Vamos agorar salvar cada uma dessas imagens para uso futuro usando:

    > library(rgdal)
    > writeRaster(I,filename='I.tif',driver='Geotiff')
    > writeRaster(H,filename='H.tif',driver='Geotiff')
    > writeRaster(S,filename='S.tif',driver='Geotiff')
    

    Parte 6 - Composição de Imagens

    Carregaremos as três imagens IHS, os primeiros cinco componentes principais e a imagem de sombra de relevo para criarmos algumas composições.

    > library(raster)
    > ihs<-stack('I.tif','H.tif','S.tif')
    > pca<-stack('pca1.tif','pca2.tif','pca3.tif','pca4.tif','pca5.tif')
    > relevo<-raster('hillshade.tif')
    

    Criaremos 4 combinações distintas com essas bandas. Experimente novas combinações.

    > par(mfrow=c(2,2))
    > PC1PC2H<-stack(pca[[1]],pca[[2]],ihs[[2]])
    > plotRGB(PC1PC2H,stretch='lin')
    > PC1PC2S<-stack(pca[[1]],pca[[2]],ihs[[3]])
    > plotRGB(PC1PC2S,stretch='lin')
    > PC1PC5H<-stack(pca[[1]],pca[[5]],ihs[[2]])
    > plotRGB(PC1PC5H,stretch='lin')
    > PC1PC5S<-stack(pca[[1]],pca[[5]],ihs[[3]])
    > plotRGB(PC1PC5S,stretch='lin')
    

    Escolhi a primeira das imagens acima que apresentou melhor resposta espectral e associei ela com a sombra de relevo. Veja o resultado final abaixo:

    > par(mfrow=c(1,1))
    > plotRGB(PC1PC2H,stretch='lin')
    > plot(relevo,add=TRUE,legend=FALSE,axes=FALSE,col=grey(0.4:0.9),alpha=0.3,zlim=c(100,200))
    

    Parte 7 - Razão entre Bandas

    Razão entre bandas é bastante usado em geoprocessamento, principalmente para se criar índices. Vamos aqui ver alguns desses índices sendo criados usando razão entre bandas.

    NDVI ou índice de vegetação normalizado resalta a diferença de resposta de clorofila em plantas nessas duas bandas, vegetação saudável (verde) tem boa reflectância para as duas bandas e vegetação seca somente para o infravermelho. O índice , através da razão, fornece valores altos para vegetação mais saudável.

    > library(raster)
    > b4<-raster('b4lr.tif')
    > b8<-raster('b8lr.tif')
    > ndvi<-(b8-b4)/(b8+b4)
    > co<-colorRampPalette(c('gray','yellow','green'), bias = 1, space = c("rgb", "Lab"),interpolate = c("linear", "spline"), alpha = FALSE)
    > plot(ndvi,col=co(24),zlim=c(-0.4,0.7))
    

    NDMI ou índice de humidade normalizado resalta a diferença de resposta da presença de humidade em vegetação saudável nessas duas bandas. O índice , através da razão, fornece valores altos para zonas com água e de maior humidade.

    > library(raster)
    > b11<-raster('b11.tif')
    > b8<-raster('b8lr.tif')
    > ndmi<-(b8-b11)/(b8+b11)
    > co<-colorRampPalette(c('gray','yellow','green'), bias = 1, space = c("rgb", "Lab"),interpolate = c("linear", "spline"), alpha = FALSE)
    > plot(ndmi,col=co(24),zlim=c(-0.4,0.7))
    

    NDWI ou índice de água normalizado resalta a presença de água nessas duas bandas. O índice , através da razão, fornece valores altos para zonas com água .

    > library(raster)
    > b5<-raster('b5.tif')
    > b3<-raster('b3lr.tif')
    > ndwi<-(b3-b5)/(b3+b5)
    > co<-colorRampPalette(c('gray','blue'), bias = 1, space = c("rgb", "Lab"),interpolate = c("linear", "spline"), alpha = FALSE)
    > plot(ndwi,col=co(24),zlim=c(0,1))