Krigagem do índice de sítio

O índice de sítio é um atributo que representa o potencial de produtividade de determinado local, sendo influenciado por fatores geológicos, topográficos, climáticos, edáficos e bióticos, apresentando de modo geral algum nível de dependência espacial, sendo portanto possível a realização de inferências espaciais sob seu comportamento.
Neste exemplo iremos interpolar o índice de sítio por meio de krigagem ordinária. A base teórica para os procedimentos abaixo descritos pode ser consultada nas teses de Mello (2004) e Pelissari (2015).

Iniciaremos com a importação de dados de unidades amostrais com índice de sítio estimado e respectivas coordenadas geográficas.

library(readxl)
library(geoR)

## Warning: package 'geoR' was built under R version 4.2.3

## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-2 (built on 2022-08-09) is now loaded
## --------------------------------------------------------------

# Importar dados
dados <- read_excel('dados_interpolacao.xlsx')
dados

## # A tibble: 101 × 4
##    PARCELA      X       Y    IS
##      <dbl>  <dbl>   <dbl> <dbl>
##  1       1 440559 7801189  31.3
##  2       2 440846 7801282  32.0
##  3       3 440734 7801475  29.6
##  4       4 440923 7801657  29.4
##  5       5 440683 7800566  33.8
##  6       6 440709 7800769  36.9
##  7       7 440528 7800929  35.7
##  8       8 440783 7800975  31.0
##  9       9 440989 7800264  30.4
## 10      10 440872 7800395  33.9
## # ℹ 91 more rows

# Transformar dados em geodata
dados <- as.geodata(dados,coords.col = 2:3, data.col = 4)

Podemos avaliar visualmente a distribuição dos dados plotando um histograma de frequência. Esta avaliação é importante para averiguarmos a necessidade de transformação nos dados, caso estes possuam uma distribuição não normal.

# Plotar histograma
hist(dados$data,
     xlab = 'Índice de Sítio (m)', ylab = 'Frequência',
     main = 'Histograma - Índice de Sítio')

Por meio de uma análise visual, os valores de sítio parecem apresentar uma distribuição próxima da normal. Vamos confirmar se há aderência aplicando o teste de Kolmogorov-Smirnov.

# Teste de Kolmogorov-Smirnov 
ks = ks.test(dados$data, 'pnorm', mean = mean(dados$data), sd = sd(dados$data))
ks

## 
## 	Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dados$data
## D = 0.050885, p-value = 0.9562
## alternative hypothesis: two-sided

O elevado p-value aponta que não há diferença estatística entre a distribuição normal e a distribuição dos dados que estamos analisando, corroborando que poderemos proceder à krigagem sem transformação da variável de interesse.
A próxima etapa será plotar e analisar o semivariograma. Para tal, utilizaremos a função variog do geoR.
Dois argumentos são importantes neste procedimento. O primeiro max.dist se refere à distância máxima para construção do variograma. É importante que este valor seja grande o suficiente para que o variograma compreenda o alcance (range) que indica a zona de influência ou de dependência espacial entre as amostras.
O segundo argumento uvec indica o número de pontos de referência (lags) que serão distribuídos proporcionalmente até atingir a distância máxima definida para o variograma.
O procedimento para obtenção de valores para estes parâmetros é iterativo e será tema de um próximo post. Por hora, definiremos o valor de 1000 para max.dist e 10 para uvec (1000/10 = 100, que é a distância aproximada entre as unidades amostrais alocadas ).

# Plotar semivariograma
semivariograma <- variog(dados, max.dist = 1000, uvec = 10,messages=FALSE)
plot(semivariograma,
     ylab = 'Semivariância', xlab = 'Distância')

O semivariograma acima considera que a semivariância possui o mesmo comportamento em todas as direções. No entanto, é fundamental avaliar se esta premissa é verdadeira. Vamos averiguar se o fenômeno da anisotropia ocorre para este conjunto de dados utilizando a função variog4.

# Investigar anisotropia
plot(variog4(dados, max.dist = 1000, uvec = 10,messages=FALSE),
     ylab = 'Semivariância', xlab = 'Distância')

Verifica-se que as linhas que representam variogramas para diferentes direções se cruzam diversas vezes, sem haver clara distinção de tendências. Deste modo, consideramos verdadeira a premissa de isotropia. Do contrário, seria necessária uma transformação das coordenadas, de modo a estabelecer a isotropia. Agora podemos proceder ao ajuste de modelos de semivariância. Neste exemplo testaremos os modelos Esférico, Exponencial, Gaussiano e Cúbico.

# Definir modelos a serem ajustados
modelos <- c('spherical','exponential','gaussian','cubic')
# Criar lista para armazenar modelos ajustados
modelos_ajustados <- list()
# Ajustar modelos de semivariância
for(i in modelos){
modelos_ajustados[[i]] <- variofit(semivariograma, cov.model = i, messages = FALSE)
}
modelos_ajustados

## $spherical
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## parameter estimates:
##    tausq  sigmasq      phi 
##   0.6878   2.4059 334.7681 
## Practical Range with cor=0.05 for asymptotic range: 334.7681
## 
## variofit: minimised weighted sum of squares = 182.0471
## 
## $exponential
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##    tausq  sigmasq      phi 
##   1.8991   1.7099 500.6367 
## Practical Range with cor=0.05 for asymptotic range: 1499.773
## 
## variofit: minimised weighted sum of squares = 130.6787
## 
## $gaussian
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##    tausq  sigmasq      phi 
##   2.3775   1.0351 557.8393 
## Practical Range with cor=0.05 for asymptotic range: 965.5186
## 
## variofit: minimised weighted sum of squares = 136.2159
## 
## $cubic
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: cubic
## parameter estimates:
##     tausq   sigmasq       phi 
##    2.3595    1.0223 1282.4003 
## Practical Range with cor=0.05 for asymptotic range: 1282.4
## 
## variofit: minimised weighted sum of squares = 135.2704

Os quatro modelos ajustados possuem três parâmetros interpretáveis: tausq se refere ao efeito pepita (nugget), sigmasq é também conhecido como partial sill e consiste no valor do patamar menos o efeito pepita, e phi se refere ao alcance (range). Vamos plotar as curvas de semivariância geradas a partir dos modelos ajustados.

par(mfrow = c(2,2))
for(i in modelos){
  plot(semivariograma,
     ylab = 'Semivariância', xlab = 'Distância',
     main = i)
  lines.variomodel(modelos_ajustados[[i]])
}

Para avaliar a eficiência dos modelos ajustados na predição do índice de sítio em locais não amostrados, devemos proceder à validação cruzada, conforme o procedimento a seguir.

# Cria lista para armazenar resultado da validação cruzada
modelos_xvalid <- list()
# Validação cruzada para todos os modelos
for(i in modelos){
modelos_xvalid[[i]] <- xvalid(dados, model = modelos_ajustados[[i]], messages = FALSE)
}
# Formatação dos resultados
erro_padrao <- unlist(lapply(modelos_xvalid,summary))[c(11,25,39,53)]
names(erro_padrao) <- modelos
desvio_do_erro <- unlist(lapply(modelos_xvalid,summary))[c(14,28,42,56)]
names(desvio_do_erro) <- modelos
# Plotar
par(mfrow = c(1,2))
text(barplot(erro_padrao, ylim =c(-0.01,0.01), main = 'Erro padrão da validação cruzada'), y = erro_padrao+sign(erro_padrao)*0.001,labels = round(erro_padrao,5))
text(barplot(desvio_do_erro, ylim = c(0,2), main = 'Desvio padrão do erro'), y = desvio_do_erro+sign(desvio_do_erro)*0.1,labels = round(desvio_do_erro,3))

As estatísticas de validação cruzada são favoráveis ao uso do modelo cúbico, portanto, utilizaremos este para interpolar o índice de sítio para toda a área e interesse. Vamos importar os limites das áreas a serem interpoladas.

# Definir modelo a ser utilizado
modelo_selecionado <- modelos_ajustados$cubic
# Importar shape com os limites da área de interesse
library(rgdal)

## Carregando pacotes exigidos: sp

## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-4, (SVN revision 1196)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/sergio.costa/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/sergio.costa/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

talhoes <- readOGR('.', 'talhoes', verbose = FALSE)

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

# Plotar talhões e pontos amostrais
points(plot(talhoes),
       x = dados[[1]][,1],y = dados[[1]][,2],
       pch='+', col='red')

Agora devemos gerar um grid de interpolação para toda a área de interesse e, por fim, aplicar o modelo cúbico ajustado.

# Gerar grid de interpolação
library(sp)
grid <- spsample(talhoes, type = 'regular', n = 10000)
plot(grid)

# Realizar a krigagem
krigagem_is <- krige.conv(dados, locations = as.data.frame(grid), krige = krige.control(obj.model = modelo_selecionado))
# Armazenar resultado em spatial data frame
krigagem_is <- coordinates(data.frame(x = data.frame(grid)[,1],
                                      y = data.frame(grid)[,2],
                                      IS = krigagem_is$predict))
# Converter em raster
library(raster)
krigagem_is <- rasterFromXYZ(krigagem_is)
# Definir projeção
projection(krigagem_is) <- CRS("+init=epsg:31982")
image(krigagem_is, asp=1,xlab='Longitude',ylab = 'Latitude')

# Plotar resultado com bordas e pontos amostrais
image(krigagem_is, asp=1,xlab='Longitude',ylab = 'Latitude')
plot(talhoes, add = TRUE)
points(x = dados[[1]][,1],y = dados[[1]][,2],
       pch='+', col='black')

Agora que visualizamos o resultado da krigagem, podemos exportar um arquivo em formato .tif, que permitirá a visualização em qualquer outro software de SIG.

# Salvar arquivo .tif
writeRaster(krigagem_is,'krigagem_is.tif', overwrite = TRUE)