Capítulo 08: Correção Geométrica e Registro de Imagens
Maximiliano Gobbo
2024-01-12
Ref: Cap. VI, itens 6.2.3.1 e 6.2.3.2 - MENESES, P. R.; ALMEIDA, T. D.; ROSA, A. N. D. C. S.; SANO, E. E. et al. Introdução ao Processamento de Imagens de Sensoriamento Remoto. Brasília: UnB-CNPq, 276 p., 2012.
1 Introdução
A correção geométrica de imagens orbitais é essencial para garantir a precisão espacial dos dados em Processamento Digital de Imagens (PDI). Este processo corrige distorções causadas pela topografia, movimento orbital e características do sensor, com os seguintes objetivos principais:
- Correção de distorções geométricas: Ajusta deformações devido à topografia e ao ângulo de visão do sensor.
- Registro espacial: Assegura que os pixels correspondam corretamente às localizações geográficas.
- Compensação de movimentos orbitais e rotação da Terra: Corrige distorções oriundas do movimento do satélite e da rotação terrestre.
- Eliminação de efeitos atmosféricos: Minimiza distorções atmosféricas na imagem.
A correção geométrica utiliza pontos de controle com coordenadas conhecidas para ajustar a imagem. Métodos comuns incluem o uso de Modelos Digitais de Elevação (MDE) e dados orbitais. No ambiente R, pacotes como raster e rgdal são utilizados para ler a imagem, identificar pontos de controle e realizar a correção geométrica. A precisão da correção depende da qualidade dos pontos de controle e do método empregado.
O processo envolve a leitura da imagem, identificação dos pontos de controle, ajuste geométrico usando pacotes específicos e reprojeção da imagem. Pontos de controle são escolhidos por sua representatividade e precisão geográfica, obtidos através de medições de campo, mapas ou imagens georreferenciadas. Ferramentas como sp, sf, raster e rgdal no R são usadas para manipular e corrigir os dados espaciais.
A eficácia da correção geométrica é avaliada pelo erro RMS, que mede a diferença entre as coordenadas estimadas e reais dos pontos de controle. O registro de imagens garante o alinhamento espacial correto, essencial para aplicações de sensoriamento remoto e análise espacial. Usar informações vetoriais, como shapefiles de limites geográficos, ajuda a orientar a correção geométrica, garantindo precisão no alinhamento da imagem com a estrutura geográfica existente.
2 Correção Geométrica e Registro de Imagem em R
Na aula de hoje iremos precisar dos pacotes a seguir. Então, antes de iniciar, realizem a instalação dos pocotes e o carregamento dos mesmos, como instruído no código a seguir.
Além disso, vamos conferir se estamos trabalhando no diretório correto (pasta de trabalho). Caso não esteja no diretório de trabalho que estamos usando, faça a alteração executando setwd().
## [1] "C:/ARQUIVOS COMPUTADOR/DOUTORADO/APOSTILA PROCESSAMENTO DIGITAL DE IMAGENS"Hoje vamos utilizar nossa imagem orbital CBERS_4_MUX_20231228_156_123_L4_BAND8.tif, mas reduziremos a resolução desta imagem para otimizar nosso processamento (tempo de máquina). O código a seguir ilustra o processo.
# Carregar a imagem orbital
imagem_orbital <- raster("CBERS_4_MUX_20231228_156_123_L4_BAND8.tif")
# Reduzir a resolução da imagem
imagem_reduzida <- aggregate(imagem_orbital, fact = 5, fun = mean) # Fator 5 para redução
# Plotar a imagem em escala de cinza
plot(imagem_reduzida, col = gray.colors(256))
Após o carregamento, vamos conferir se a imagem carregamento passou por algum pré-processamento, ou seja, se foi inserido nela algum sistema de coordenadas.
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["UTM zone 23S",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-45,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",10000000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",17023]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]Com base no resultado, observamos que a imagem já possui um sistema de coordenadas associado. O sistema de coordenadas é presentado em UTM, zona 23S, com o datum WGS84. Diante disso, vamos realizar dois ajustes geométricos. O primeiro será a substituição do sistema de coordenadas pelos SIRGAS2000, sem nenhum ponto de correção. No segundo, realizaremos a correção geométrica utilizando o shapefile dos municípios paulistas, tendo em vista que a imagem orbital focaliza-se nesta região.
Para a correção geométrica através da substituição do sistema de coordenadas, seguiremos o código a seguir.
# Definir o sistema de coordenadas SIRGAS 2000 usando o código EPSG
novo_crs <- CRS("+init=EPSG:31983")
# Projetar a imagem para o novo sistema de coordenadas
imagem_sirgas2000 <- projectRaster(imagem_reduzida, crs = novo_crs)
# Verificar o sistema de coordenadas da imagem corrigida
crs(imagem_sirgas2000)## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["SIRGAS 2000 / UTM zone 23S",
## BASEGEOGCRS["SIRGAS 2000",
## DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4674]],
## CONVERSION["UTM zone 23S",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-45,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",10000000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",16123]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## USAGE[
## SCOPE["unknown"],
## AREA["Brazil - between 48°W and 42°W, northern and southern hemispheres, onshore and offshore."],
## BBOX[-33.5,-48,5.13,-42]]]O código acima mostra que o sistema de coordenadas SIRGAS2000 foi aplicado à imagem orbital, portanto, nossa imagem_sirgas2000 foi reprojetada no plano segundo este sistema de referenciamento.Vamos comparar visualmente se a alteração do sistema de coordenadas da imagem refletiu significativamente na projeção desta imagem. O código a seguir realiza este procedimento.
par(mfrow=c(1,2))
plot(imagem_reduzida, col = gray.colors(256), main = "Imagem Reduzida")
plot(imagem_sirgas2000, col = gray.colors(256), main = "Imagem Corrigida (SIRGAS 2000)")
Conforme observamos, não houve mudança significativa na projeção da plot_sirgas2000 quando comparada a plot_reduzida. Assim, vamos prosseguir a com a correção geométrica e adentrar nosso segundo passo.
Neste momento precisamos estabelecer os pontos de correção, os quais representam pontos conhecidos na superfície para os quais conhecemos as coordenadas geográficas reais e que são empregados no ajuste geométrico da imagem orbital, de modo a alinhar a imagem orbital corretamente com as coordenadas geográficas.
Esses pontos podem ser obtidos de diferentes fontes, como mapas georreferenciados, dados de GPS ou pontos de controle coletados em campo. Uma vez que se tenha os pontos de correção, pode-se usá-los juntamente com as funções em pacotes rgdal para realizar a correção geométrica.
O shapefile pode ser usado como referência para selecionar e distribuir os pontos de correção na imagem orbital. Por exemplo, se o shapefile contiver informações sobre limites administrativos, como fronteiras municipais ou estaduais, você pode selecionar pontos de correção em locais estratégicos dentro desses limites.
Além disso, se você tiver coordenadas geográficas conhecidas para esses pontos no shapefile, poderá usá-las diretamente como pontos de correção. Isso ajudará a garantir que a imagem orbital seja corrigida e alinhada corretamente com a área de interesse representada pelo shapefile.
Portanto, o shapefile pode servir como uma ferramenta valiosa para orientar o processo de correção geométrica, fornecendo pontos de referência geográficos que podem ser usados para alinhar adequadamente a imagem orbital.
# Definir o caminho para o arquivo shapefile (.shp)
cidades_paulistas <- "sp_municipios.shp"
# Carregar o shapefile
dados_cidades_paulistas <- st_read(cidades_paulistas)## Reading layer `sp_municipios' from data source
## `C:\ARQUIVOS COMPUTADOR\DOUTORADO\APOSTILA PROCESSAMENTO DIGITAL DE IMAGENS\sp_municipios.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 645 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -53.11011 ymin: -25.31232 xmax: -44.16137 ymax: -19.77966
## Geodetic CRS: SIRGAS 2000## Simple feature collection with 645 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -53.11011 ymin: -25.31232 xmax: -44.16137 ymax: -19.77966
## Geodetic CRS: SIRGAS 2000
## First 10 features:
## ID CD_GEOCODM NM_MUNICIP geometry
## 1 1727 3500105 ADAMANTINA MULTIPOLYGON (((-51.09093 -...
## 2 1728 3500204 ADOLFO MULTIPOLYGON (((-49.69668 -...
## 3 1729 3500303 AGUA\xcd MULTIPOLYGON (((-47.01254 -...
## 4 1730 3500402 \xc1GUAS DA PRATA MULTIPOLYGON (((-46.73069 -...
## 5 1731 3500501 \xc1GUAS DE LIND\xd3IA MULTIPOLYGON (((-46.635 -22...
## 6 1732 3500550 \xc1GUAS DE SANTA B\xc1RBARA MULTIPOLYGON (((-49.28903 -...
## 7 1733 3500600 \xc1GUAS DE S\xc3O PEDRO MULTIPOLYGON (((-47.86292 -...
## 8 1734 3500709 AGUDOS MULTIPOLYGON (((-49.0167 -2...
## 9 1735 3500758 ALAMBARI MULTIPOLYGON (((-47.86567 -...
## 10 1736 3500808 ALFREDO MARCONDES MULTIPOLYGON (((-51.37102 -...
Com os dados shapefile referente aos municípios paulistas carregados, vamos realizar uma filtragem por municípios que estão na nossa área de interesse, o que chamaremos de entorno_franca, tendo em vista que a cidade de Franca encontra-se no centro da imagem orbital. Acompanhe no código a seguir.
# Filtrar o shapefile para selecionar apenas os dados do entorno de Franca
entorno_franca <- subset(dados_cidades_paulistas, NM_MUNICIP %in% c("BATATAIS", "BURITIZAL", "CRISTAIS PAULISTA", "FRANCA", "IGARAPAVA", "ITUVERAVA", "JERIQUARA", "PEDREGULHO", "RESTINGA", "RIFAINA"))
# Visualizar os dados da cidade de Franca
print(entorno_franca)## Simple feature collection with 10 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -48.04027 ymin: -21.03843 xmax: -47.14346 ymax: -19.9607
## Geodetic CRS: SIRGAS 2000
## ID CD_GEOCODM NM_MUNICIP geometry
## 67 1793 3505906 BATATAIS MULTIPOLYGON (((-47.43263 -...
## 94 1820 3508207 BURITIZAL MULTIPOLYGON (((-47.64544 -...
## 147 1873 3513207 CRISTAIS PAULISTA MULTIPOLYGON (((-47.36814 -...
## 186 1912 3516200 FRANCA MULTIPOLYGON (((-47.41682 -...
## 230 1956 3520103 IGARAPAVA MULTIPOLYGON (((-47.86035 -...
## 275 2001 3524105 ITUVERAVA MULTIPOLYGON (((-47.9233 -2...
## 288 2014 3525409 JERIQUARA MULTIPOLYGON (((-47.62619 -...
## 417 2143 3537008 PEDREGULHO MULTIPOLYGON (((-47.57243 -...
## 479 2205 3542701 RESTINGA MULTIPOLYGON (((-47.44607 -...
## 490 2216 3543600 RIFAINA MULTIPOLYGON (((-47.35998 -...Para selecionar os pontos de correção, podemos usar os vértices das geometrias das cidades próximas a Franca. Para extrair os vértices das geometrias das cidades, podemos usar a função st_coordinates do pacote sf. Esta função extrai as coordenadas dos vértices de uma geometria. Posteriormente, precisamos convertê-los para as coordenadas do sistema de referência da imagem orbital antes de aplicá-los como pontos de controle.
# Extrair os vértices das geometrias das cidades próximas a Franca
vertices_entorno <- st_coordinates(entorno_franca$geometry)
# Criar um objeto sf diretamente a partir dos vértices #Defina o CRS corretamente conforme necessário
vertices_entorno_sf <- st_as_sf(data.frame(lon = vertices_entorno[,1], lat = vertices_entorno[,2]), coords = c("lon", "lat"), crs = st_crs("EPSG:31983")) # Extrair os vértices do raster imagem_sirgas2000
vertices_imagemsirgas <- rasterToPoints(imagem_sirgas2000)
# Criar um data frame com os vértices
vertices_df <- data.frame(X = vertices_imagemsirgas[, 1], Y = vertices_imagemsirgas[, 2])
# Converter para um objeto sf
vertices_imagemsirgas_sf <- st_as_sf(vertices_df, coords = c("X", "Y"), crs = st_crs(imagem_sirgas2000))Agora que temos os vértices convertidos em um objeto sf, podemos prosseguir para verificar quais desses pontos de vértice estão contidos na extensão da imagem orbital. Mas, antes de avançarmos, vamos comparar se vertices_entorno_sf e vertices_imagemsirgas_sf possuem o mesmo sistema de coordenadas geográficas e se possuem a mesma extensão (XY mínimo e máximo).
# Transformar o sistema de coordenadas dos pontos de controle para o mesmo sistema
vertices_entorno_sf_trans <- st_transform(vertices_entorno_sf, crs = st_crs(vertices_imagemsirgas_sf))
# Comparar CRS
if (identical(st_crs(vertices_entorno_sf_trans), st_crs(vertices_imagemsirgas_sf))) {
cat("Ambos os dataframes possuem o mesmo CRS.\n")
} else {
cat("Os dataframes possuem CRS diferentes.\n")
}## Ambos os dataframes possuem o mesmo CRS.Conferido os crs dos conjuntos passamos à junção dos pontos correspondentes entre a nossa referência (vertices_entorno_sf_trans) e o nosso alvo (vertices_imagemsirgas_sf). Na sequência, iremos extrair as coordenadas do nosso conjunto pontos_correspondentes, que são os dados que possuem mesma posição na imagem.
# Encontrar os pontos correspondentes nos conjuntos de dados de referência e alvo
pontos_correspondentes <- st_join(vertices_entorno_sf_trans, vertices_imagemsirgas_sf, join = st_nearest_feature)Extraídas as coordenadas, passaremos à “construção do modelo linear” (2D) para correção da nossa imagem. O código a seguir ilustra o procedimento de construção (ajuste) do modelo linear, a extração dos coeficientes do modelo e a construção da matriz de transformação para correção geométrica. Após a construção da matriz de transformação, definiremos o crs para a imagem corrigida e aplicaremos a matriz de transformação para corrigir nossa imagem orbital.
# Ajustar um modelo linear aos pontos de controle
modelo <- lm(coordenadas[, "Y"] ~ coordenadas[, "X"])
# Extrair os coeficientes do modelo
coeficientes <- coef(modelo)
# Extrair os coeficientes do modelo
a <- coeficientes[1]
b <- coeficientes[2]
# Construir a matriz de transformação
transformacao <- matrix(c(a, b, 0, a), nrow = 2, byrow = TRUE)
# Definir a matriz de transformação
matriz_afinidade <- matrix(c(coeficientes[2], coeficientes[4], 0, coeficientes[1], coeficientes[3], 0), ncol = 2, byrow = TRUE)
# Definir o CRS da imagem corrigida
crs_corrigida <- CRS("+proj=utm +zone=23 +datum=WGS84 +units=m +no_defs")
# Aplicar a transformação à imagem alvo
imagem_corrigida <- projectRaster(imagem_sirgas2000, crs = crs_corrigida, transform = matriz_afinidade)Após a aplicação da matriz de transformação, vamos analisar visualmente como ficaram nossas imagens, antes e depois da correção geométrica.
par(mfrow=c(1,2))
plot(imagem_corrigida, col = gray.colors(256), main = "Imagem com Correcao Geometrica")
plot(imagem_sirgas2000, col = gray.colors(256), main = "Imagem SIRGAS 2000")
No contexto de processamento de imagens, o registro de imagens refere-se ao processo de alinhar duas ou mais imagens para que correspondam geometricamente umas às outras. Isso é importante em muitas aplicações, como em análise de mudanças ao longo do tempo, fusão de imagens de diferentes sensores ou fontes, e sobreposição de dados de diferentes modalidades para análises integradas.
No nosso caso, você corrigimos geometricamente uma imagem para se alinhar com uma imagem de referência (SIRGAS 2000). Esse processo de correção geométrica pode ser considerado um tipo de registro de imagem, onde a imagem corrigida foi ajustada para se sobrepor adequadamente à imagem de referência. Isso é importante para garantir que as duas imagens estejam alinhadas corretamente, o que facilita a comparação e a análise conjunta das informações contidas nelas. E assim encerramos nossa prática sobre correção geométrica e registro de imagem.
Capítulo 07 - Aritmética de Bandas: soma, subtração, multiplicação e divisão
3 Referências Bibliográficas
ANDERSON, E. C. Making Maps with R Â: Reproducible Research. GitHub, 2021. Disponível em: https://eriqande.github.io/rep-res-web/lectures/making-maps-with-R.html.
CHAMBERS, J. M. Extending R. CRC Press, 364p., 2017. ISBN: 9781498775724.
DA SILVA, F. R.; GONÇALVES-SOUZA, T.; PATERNO, G. B.; PROVETE, D. B.; VANCINE, M. H. Análises ecológicas no R. Recife, PE: Nupeea, 1ª ed., 640p., 2022. ISBN: 9788579175640. Disponível em: https://analises-ecologicas.com/.
GAUTHIER, N. Anthromes 12K DGG (V1) analysis code and R research compendium. Harvard Dataverse, 2021. DOI: https://doi.org/10.7910/DVN/6FWPZ9.
LOVELACE, R. Introduction to visualizing spatial data in R. 2017. Disponível em:https://github.com/Robinlovelace/Creating-maps-in-R.
LOVELACE, R., NOWOSAD, J., & MUENCHOW, J. Geocomputation with R. Chapman and Hall/CRC, 353p., 2019. Disponível em: https://www.routledge.com/Geocomputation-withR/Lovelace-Nowosad-Muenchow/p/book/9780367670573
MACHLIS, S. Create maps in R in 10 (fairly) easy steps. Computerworld, 2017. Disponível em: https://www.computerworld.com/article/3038270/create-maps-in-r-in-10-fairly-easy-steps.html.
MARTOS, G. Cluster Analysis with R. 2021. Disponível em: https://rstudio-pubsstatic.s3.amazonaws.com/33876_1d7794d9a86647ca90c4f182df93f0e8.html.
MENESES, P. R.; ALMEIDA, T. D.; ROSA, A. N. D. C. S.; SANO, E. E. et al. Introdução ao Processamento de Imagens de Sensoriamento Remoto. Brasília: UnB-CNPq, 276 p., 2012.
MORENO, M. & BASILLE, M. Drawing beautiful maps programmatically with R, sf and ggplot2 - Part 1: Basics. 2021. Disponível em: https://r-spatial.org/r/2018/10/25/ggplot2-sf.html.
MUENCHOW, J.; SCHRATZ, P.; BRENNING, A. RQGIS: Integrating R with QGIS for Statistical Geocomputing. The R Journal, 9, n. 2, p. 409-428, 2017 2017. Doi: https://doi.org/10.32614/RJ-2017-067.
NOVIA, D. Cluster Analysis in R: Tips for Great Analysis and Visualization. Datanovia, 2018.Disponível em: https://www.datanovia.com/en/blog/cluster-analysis-in-r-simplified-and-enhanced/.
PEBESMA, E.; BIVAND, R. Spatial Data Science: With Applications in R. Chapman and Hall/CRC, 1ª ed., 314p., 2023. Doi: https://doi.org/10.1201/9780429459016.
R-BLOGGERS. Cluster Analysis in R. 2021. Disponível em: https://www.rbloggers.com/2021/04/cluster-analysis-in-r/.
STANIAK, M.; BIECEK, P. The Landscape of R Packages for Automated Exploratory Data Analysis. The R Journal, v. 11, n. 2, 23 p., 2019. Disponível em: https://arxiv.org/pdf/1904.02101.
VENABLES, W. N.; SMITH, D. M. An Introduction to R. R Core Team, 2017. Disponível em: http://colinfay.me/intro-to-r/#.
WICKHAM, H. Advanced R. CRC Press, 2ª ed., 2022. Disponível em: https://advr.hadley.nz/#license
:::