Cálculo de Shannon ' diversidad s ventana móvil de R

Question

Tengo un mapa de bits, donde cada pixel es un determinado tipo de vegetación. Estoy interesado en la diversidad de tipos de vegetación alrededor de cada píxel.

Me gustaría asignar un índice de diversidad de Shannon (ver este enlace para la ecuación) para cada píxel en función de una ventana móvil de análisis de los píxeles circundantes.

He llegado con un poco de código para contar el número de valores únicos en el movimiento de la ventana de la utilización de funciones en la raster paquete (al menos eso es lo que espero que se está haciendo!), pero estoy teniendo problemas con el resto del código necesario. Aquí es lo que tengo hasta ahora:

fw<-focalWeight(veg_map,5000, "circle") # creates circular filter with a radius of 5000m
test_fcl<-focal(veg_map,w=fw, fun=function(x,
          ...){length(unique(na.omit(x))) }) #counts unique values in moving window #  (e.g. species richness)

La parte estoy teniendo más problemas con es encontrar la proporción de cada tipo único dentro de la ventana móvil. Desde que yo podría calcular varias métricas de diversidad, incluyendo la de Shannon.

Answer 1

Calcular focal medio de los indicadores de cada tipo de vegetación. En cada celda, estos dan las proporciones de los tipos. Multiplicar cada uno por su logaritmo negativo y suma: es el índice de diversidad.

Usted encontrará que incluso para un gran número de categorías (incluso cientos), esto es increíblemente más rápido que el método de fuerza bruta de las tablas de cada barrio en turno, incluso con pequeños barrios. Esto es debido a que se lleva a cabo por medio de las circunvoluciones, el logro de O(m*n*log(m*n)) ajuste de la escala de m filas y n columnas, independientemente de la vecindad de tamaño. La fuerza bruta método de tabulación simplemente no va a completar la ejecución en cualquier momento razonable, en cualquier rásteres de tamaño significativo o el uso de los grandes barrios.

El código al final se ilustra el procedimiento. Se utiliza un 41 41 circular barrio en un 500 por 800 trama que incluye cinco categorías; el tiempo es menos de 10 segundos. Que es pobre, en realidad: puede ser mejorado por varios órdenes de magnitud en otras plataformas. Sin embargo, podría ser lo suficientemente bueno para el trabajo de producción cuando los rásteres no son demasiado grandes.

Este código da resultados que difieren de vegan y de otra respuesta en este hilo. La razón es sutil: focal primer multiplica todos los elementos en el barrio por el centro de pesos y, a continuación, pasa todos los valores a su función. Para no rectangular barrios que presentan algunos de los ceros. Los ceros no están incluidas en el barrio y por lo tanto no deben ser tabulados. Sin embargo, el enfoque utilizando table realmente hace el cómputo de las frecuencias de los ceros. El efecto es añadir un valor constante (igual a r *log(r) donde r es la proporción de ceros en un barrio), excepto alrededor de los bordes (donde el barrio de la forma, y por lo tanto r, puede cambiar).

library(raster)
m <- 500            # Rows
n <- 800            # Columns
cellsize <- 250     # Meters
p <- (1:5)^2        # Relative probabilities of vegetation; first is NoData
radius <- 5000      # Meters
#
# Create sample data.
#
set.seed(17)
x <- matrix(sample.int(length(p), m*n, replace=TRUE, prob=p), nrow=m)
#
# Convert to raster.
#
x.r <- raster(x, xmn=0, ymn=0, xmx=n*cellsize, ymx=m*cellsize)
#
# Diversity index.
#
diversity <- function(x.r, radius, type="circle", verbose=TRUE) {
  # 
  # Create focal weights matrix.
  #
  nbrhd <- focalWeight(x.r, radius, type=type)
  #
  # Compute focal means of indicators.
  #
  entropy.r <- calc(x.r, function(x) 0)
  x.log <- function(x) ifelse(x==0, 0, x*log(x))
  for (i in 1:length(p)) {
    if (verbose) cat("Computing indicator for category", i, "... ")
    z <- system.time({
      entropy.r <- entropy.r - calc(focal(x.r == i, nbrhd), x.log)
    })
    if (verbose) cat(z[3], "seconds.\n")
  }
  return (entropy.r)
}
#
# Plot the original and its entropy.
#
entropy.r <- diversity(x.r, radius)
par(mfrow=c(1,2))
plot(x.r, main="Original")
plot(entropy.r, main="Entropy")

Answer 2

He aquí una sencilla aplicación para usted. Edit: focal fija de pesos de la matriz a excluir 0s como por whuber comentarios en su respuesta.

library(raster)

# Example Data
set.seed(1)
r <- raster(matrix(sample(1:10, 100, replace=T), 10, 10))

# Calculate a weights matrix, and reset elements to 0s and 1s 
# rather than true weights
fw <- focalWeight(r, 0.2, 'circle')
fw <- ifelse(fw == 0, NA, 1)

# Neighbourhood richness
richness <- function(x, ...) {
    length(unique(na.omit(x)))
}
richOut <- focal(r, fw, fun=richness, pad=T)

# Neighbourhood Shannon Diversity Index
shannon <- function(x, ...) {
    cnts <- table(x)
    cnts <- cnts / sum(cnts)
    -sum(cnts * log(cnts))
}
shanOut <- focal(r, fw, fun=shannon, pad=T)

El vegan paquete también tiene un diversity función que puede calcular Shannon, Simpson, y los índices de diversidad de Fisher para usted. Los resultados son los mismos.

library(vegan)
shannonVegan <- function(x, ...) {
    diversity(table(x), index="shannon")
}
shanVegOut <- focal(r, fw, fun=shannon, pad=T)  

Se me ocurre que puede haber algunos efectos de borde debido a la focal círculo de ir fuera del borde de la cuadrícula. Usted podría cambiar pad=T a pad=F conseguir NAs para ninguna de las celdas donde esto sería un problema, o simplemente ser consciente de ello.