Create presence/absence (0/1) rasters for every value in categorical raster in R?

Question

I am conducting a logistic regression analysis of spatial data for a project using STAN MCMC within R. The regression analysis itself is done, but now I need a way to apply the parameter estimates for different variables to raster data to generate prediction maps. Applying parameter estimates to continuous rasters in an ArcGIS Raster Calculator operation is straightforward (e.g., Intercept + 0.44*elevation - 0.23*dist_to_water...), but incorporating categorical rasters is more difficult. STAN requires categorical variables to be split up into a series of dummy variables, so my categorical rasters (e.g., native veg, surface geology, erosion class) need to be split up into a series of presence/absence (0/1) rasters for each value. For example, Native Vegetation has four categories, so I need four different 0/1 rasters that correspond to the distribution of each vegetation type. This is a similar problem to that discussed here.

However, I have six categorical rasters that each include a large number of values. As such, I would like to avoid having to manually create 0/1 rasters for each value using tools in ArcGIS like Extract by Attributes, Reclassify, or Conditional statements (e.g., Con(native_veg, 1, 0, "VALUE = 1")).

Is there a function in R that can iterate through unique raster values and create new 0/1 rasters for each?

Answer 1

The raster package has a function to do this for you in one line. Use layerize():

library(raster)
# make an example raster
r <- raster(nrow=100, ncol=100)
r[] <- round(runif(ncell(r),1,4),0)

# create presence/absence rasters, stored in a RasterBrick. 
r_dummy <- layerize(r)
plot(r_dummy[[1]])

Answer 2

Are you sure that the dummy coding in STAN requires separate binomial parameters for the estimate? For efficiency sake, the dummy coding for say, Native Vegetation should be a factor of [1,2,3,4] and not four separate parameters. If this is, in fact, the case it is trivial to do this in R.

library(raster)
r <- raster(nrow=100, ncol=100)
  r[] <- round(runif(ncell(r),1,4),0)

r1=r; r2=r; r3=r; r4=r # make copies of original to modify
  r1[] <- ifelse(r[] == 1, 1, 0)
  r2[] <- ifelse(r[] == 2, 1, 0)
  r3[] <- ifelse(r[] == 3, 1, 0)
  r4[] <- ifelse(r[] == 4, 1, 0)

plot(stack(r1,r2,r3,r4))

You can extend this into processing all of your rasters in a for loop.

First, as an example, create raster stack with 5 layers and different unique values ranging 1 - 5. This should emulate you problem since each raster has an different number of levels. This simulated data is, of course, representing a raster stack or brick object of you nominal covariates that need to be converted to binomial. Your raster data can be read into this object class using the raster::stack or raster::brick functions.

r <- stack(brick(array(runif(100 * 100 * 5), dim=c(100, 100, 5))))
  for(i in 1:nlayers(r)) { r[[i]][] <- round(runif(ncell(r), 
                           sample(1:2,1),sample(3:5,1) ),0) }
names(r) <- paste0("parameter",1:5)

Now we can define a double loop that looks at the unique values for each raster and then loops through them to create binary rasters for each rasters unique levels. The resulting object (given the below code) will be binary.rasters.

binary.rasters <- stack()
rnames <- vector() 
  for(i in 1:nlayers(r)) {
    for(j in 1:length(unique(r[[i]][]))) {
      u <- unique(r[[i]][])
      b <- r[[i]] 
      b[] <- ifelse(b[] == u[j], 1, 0)
      binary.rasters <- addLayer(binary.rasters, b)
      rnames <- append(rnames, paste(names(r)[i], paste("level",u[j],sep="-"), sep="_") )  
    }
  }
( names(binary.rasters) <- rnames )

I added a vector rnames that tracks the raster name and appends it with the level in each sub loop. This can then be added as the names of the resulting raster stack. This will allow one to know what raster and level a given binary raster is representing however, does not necessary correspond with the parameter names used in the original model. Although, you can easily rename the elements in this vector.

I have seen this type of dummy coding in occupancy modeling (eg., software Presence) and it sure eats up parameter space. One has to wonder if it would not be prudent to form explicit hypothesis around a specific level in a given categorical variable and not think of it as a single independent variable, because that is sure not how it is behaving in the model. In this case each level in the variable is truly a separate parameter. If there are only a level or two in a given independent nominal variable that you hypothesize is going to effect your process, why include all of the levels? By reducing the parameter space you are going to get much more relevant AIC values when evaluating competing models.