Sen Slope para Sensoriamento Remoto

Antonio Barros

O estimador de Theil-Sen de um conjunto de pontos bidimensionais (\(x_i\), \(y_i\)) é a mediana m das inclinações \({(y_j – y_i)}/{(x_j – x_i)}\) determinadas por todos os pares de pontos amostrais. Sen (1968) estendeu essa definição para lidar com o caso em que dois pontos de dados têm a mesma coordenada x. Na definição de Sen, toma-se a mediana das inclinações definidas apenas a partir de pares de pontos com coordenadas x distintas.

\[ \text{Sen's Slope} = \text{mediana}(x_j - x_i), \forall i < j \] em que \(x_j - x_i\) são os dados na série histórica do tempo i e j.

Considerando n número de pontos dos dados, então

\[N = \frac{n(n-1)}{2}\] onde N é o o número calculado de estimativas da inclinação.

A inclinação de Sen’s slope é dada como a mediana do número de N por meio da inclinação. Se N é par \([(N+1)/2]\), caso seja ímpar \([(N/2) + ((n/2)/2)]/2\)

SenSlopeRaster <- function(rasterstack) {
  print(paste("Iniciando SenSlopeRaster:", Sys.time()))
  layers <- nlayers(rasterstack)
  ncell <- ncell(rasterstack)
  ncol <- ncol(rasterstack)
  nrow <- nrow(rasterstack)
  crs <- crs(rasterstack)
  extent <- extent(rasterstack)
  pb <- txtProgressBar(min = 0, max = ncell, initial = 0)
  print("Parâmetros carregados")
  
  mtrx <- as.matrix(rasterstack, ncol = layers)
  empt <- matrix(nrow = ncell, ncol = 1)
  print("Iniciando loop de cálculo do Sen Slope")
  for (i in 1:ncell) {
    y <- mtrx[i, ]
    
    if (all(is.na(y)) || sum(!is.na(y)) < 3) {
      empt[i, 1] <- NA
    } else {
      y <- na.omit(y)
      SenSlope <- sens.slope(y)$estimates
      empt[i, 1] <- SenSlope
    }
    setTxtProgressBar(pb, i)
  }
  
  SenSlope_raster <- raster(nrows = nrow, ncols = ncol, crs = crs)
  extent(SenSlope_raster) <- extent
  values(SenSlope_raster) <- empt[, 1]

  print(paste("Sen Slope calculado com sucesso. Fim da função em", Sys.time()))
  return(SenSlope_raster)
}

Como executar:

#sen_slope <- SenSlopeRaster(pilha_raster)

Testando:

## Carregando pacotes exigidos: sp

# Calculando o Sen Slope para cada pixel usando a função SenSlopeRaster
sen_slope <- SenSlopeRaster(Conjunto)

## [1] "Iniciando SenSlopeRaster: 2023-05-24 18:14:59.070628"
## [1] "Parâmetros carregados"
## [1] "Iniciando loop de cálculo do Sen Slope"
## ================================================================================[1] "Sen Slope calculado com sucesso. Fim da função em 2023-05-24 18:15:54.97461"

# Resultado
plot(sen_slope)