I have a raster where each pixel is a particular vegetation type. I'm interested in the diversity of vegetation types around each pixel.

I'd like to assign a Shannon's diversity index value (see this link for equation) to each pixel based on a moving window analysis of surrounding pixels.

I've come up with some code to count the number of unique values in the moving window using functions in the raster package (at least that is what I hope it is doing!), but I am having trouble with the rest of the needed code. Here is what I have so far:

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)

The part I'm having the most trouble with is finding the proportion of each unique type within the moving window. From that I could calculate several metrics of diversity, including Shannon's.

  • 1
    @aaryno. My pixels are 60 x 60m. My understanding from the Raster manual is that if the filter is circular is uses map units. Is that incorrect? Even still, just to count the unique codes takes a bit of time due to the size of my landscape. – KevinB Jun 23 '15 at 20:05
  • @KevinB - my bad - you are right, the radius is calculate in units of the CRS. to clear things up I'm deleting my comment. – aaryno Jun 23 '15 at 21:11
up vote 7 down vote accepted

Compute focal means of the indicators of each vegetation type. At each cell, these give the proportions of the types. Multiply each by its negative logarithm and sum: that's the diversity index.

You will find that even for large numbers of categories (even into the hundreds), this is fantastically faster than the brute-force method of tabulating each neighborhood in turn, even with tiny neighborhoods. That is because it is carried out by means of convolutions, achieving O(m*n*log(m*n)) scaling for m rows and n columns, regardless of the neighborhood size. The brute-force tabulation method simply will not complete executing in any reasonable time on any rasters of meaningful size or using large neighborhoods.

The code at the end illustrates the procedure. It uses a 41 by 41 circular neighborhood on a 500 by 800 raster involving five categories; the timing is under 10 seconds. That's poor, actually: it can be improved by several of orders of magnitude on other platforms. However, it might be good enough for production work when rasters are not too large.

This code gives results that differ from vegan and from another answer in this thread. The reason is subtle: focal first multiplies all elements in the neighborhood by the focal weights and then passes all values to its function. For non-rectangular neighborhoods that introduces some zeros. Those zeros are not included in the neighborhood and therefore should not be tabulated. However, the approach using table actually does tabulate the frequencies of those zeros. The effect is to add a constant value (equal to -r *log(r) where r is the proportion of zeros in a neighborhood) except around the edges (where the neighborhood shape, and therefore r, can change).

Figure

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")
  • thanks for the code. It appears to do exactly what I was hoping to do. One question, when I was first testing it I retained the original pixel values, which ranged from 0 to 20000 (though there were only 13 unique values in that range). When I ran it with those values the result would be 0 across the entire raster. I went back and changed the values to be in sequential order (i.e. 0-13) and the code seemed to function correctly. Why did the value of the pixel matter for the calculation? – KevinB Jun 26 '15 at 21:46
  • 1
    I voted this response as the answer because not only did it provide an solution but was also did so in a reasonable time (for large landscapes) compared to the other answer that was provided. – KevinB Jun 29 '15 at 16:28
  • I suspect your problem might have been that the example code assumes the pixel values are the same as the indexes in 1:length(p). When that is not the case, you need to loop over the values themselves rather than over the indexes. Nothing else needs to change. – whuber Jun 29 '15 at 17:35
  • I tried to run this over an NDVI image, but always get a raster with only one value or a small spot on a small part of the raster. Is there any way to make this work over a continuous raster dataset? – mace May 27 at 7:18
  • @mace That reflects a problem in your code rather than a problem with this algorithm. Indicator functions don't care whether values are "continuous" or discrete or colors or anything else. – whuber May 27 at 11:38

Here's a simple implementation for you. Edit: fixed focal weights matrix to exclude 0s as per whuber's comments in his answer.

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)

The vegan package also has a diversity function that can calculate Shannon, Simpson, and Fisher diversity indices for you. The results are the same.

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

It occurs to me that there may be some edge effects due to the focal circle going off the edge of the grid. You could change pad=T to pad=F to get NAs for any cells where this would be an issue, or just be aware of it.
Imgur

  • What I love about this answer is it defines the functions for richness and diversity so clearly and understandably, making the call to focal very readable. – aaryno Jun 23 '15 at 21:06
  • @Matt SM- Thanks! I think this will work. I'll have to think about using pad. I've found that when pad =T that it creates artificial edge effects. – KevinB Jun 23 '15 at 21:25
  • 1
    I do like the clarity, @aaryno, and therefore was quite disturbed when I discovered that my code gave different answers than this solution. I delayed posting my answer until I could track down the reason (by manually computing the values for simple small grids), which I have explained in that answer. It comes down to a misunderstanding of how focal works. I have to conclude that vegan--and therefore also this answer--gives invalid values. – whuber Jun 23 '15 at 21:40
  • 1
    Thank you for showing the workaround with the focal weights. I am curious about the timing for your calculations and especially their scalability. – whuber Jun 23 '15 at 21:50
  • 1
    So, are you not biasing towards alpha diversity by calculating the index in a local window without accounting for the population? In this case I believe that "p" should represent a proportion based on the number of possible classes and not the number of classes available to the focal function. Easy fix, just modify the diversity function to account for all the classes in the raster. – Jeffrey Evans Jun 24 '15 at 6:51

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.