Generating random raster using weights with negative binomial distribution using R?

Question

I'm trying to generate "species counts" in raster form -- essentially, the number of individuals in a species in a raster cell. I'd like to use the negative binomial distribution to do so. Normally I would just use something like the Create Random Raster tool in ArcGIS Desktop; however, there is an extra wrinkle. I'd also like the species to be weighted based on a second raster, where higher values are more likely to fall in certain areas.

So if I were hypothetically modeling a species like a lion, I'd create a raster with possible values taken from a negative binomial distribution where higher values were more likely on savannahs, less likely in urban areas, and never found in the middle of lakes or oceans (with the land classes taken from an NDVI raster, or something similar).

Is the R raster package capable of doing this?

Answer 1

Using these packages (MASS gets you the negative binomial distribution functions):

> library(MASS)
> library(raster)

Let's create a raster of coefficients that you have computed in some way from your landscape. Here its a 3x4 grid of numbers from 1 to 12:

> C = raster(matrix(1:12,3,4))

Then we can create negative binomial samples with a mu value matching the values in C.

First create an empty raster of the same size and shape:

> P = raster(matrix(NA, 3, 4))

Then replace the values using rnegbin with the mu vector from the corresponding values in C. Here I've set theta to a constant 1.

> P[] = rnegbin(ncell(C), mu=C[], theta=1)
> P[]
 [1]  0  0  0  8  0  1 10  3  3  5  2  8

Plotting P will show its a very noisy version of C - there's a lot of zeroes and you only start getting positive numbers for larger values of C. Such is randomness.

Depending on how you want to parameterise your negative binomial distributed samples you might want to set theta from another raster, or use the parameterisation of the negative binomial from the rnbinom function.