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)
## Warning: package 'readxl' was built under R version 3.6.1
library(geoR)
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------
# Importar dados
dados <- read_excel('dados.xlsx')
dados
## # A tibble: 101 x 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
## # ... with 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.

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

O gráfico indica que os valores de sítio apresentam uma distribuição próxima da normal, o que significa que não há necessidade de transformação dos dados. Vamos confirmar 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
## 
##  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. Realizadas as devidas análises, procedemos 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)
## Loading required package: sp
## rgdal: version: 1.4-4, (SVN revision 833)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/sergio.costa/Documents/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/sergio.costa/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.3-1
talhoes <- readOGR('.', 'talhoes', verbose = FALSE)
# 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)

Espero que este post lhe tenha sido útil!
Qualquer dúvida, entre em contato pelo e-mail: