Cross tab of two rasters across different polygons in a shapefile

Question

I need to find the number of pixels of raster 1 that fall within different pixel values of raster 2. For example, there are 250 pixels of value=2 in raster1 that in value=89 of raster 2. Hence, I am using crosstab() from the raster package.

Now I need to find the cross tab of pixel count from the two rasters, within polygons in a shapefile. How do I do that?

One can think of this problem as- area (pixel count* pixel size) of opportunity (0-no opportunity, 1-opportunity specified in pixels in raster 1) in different categories of standards of living (1- low standard of living, 2-medium, 3-high specified in pixels in raster 2), across different US counties/states (polygons in shapefile).

Answer 1

Rasterize the polygons to the same basis as your rasters, then use crosstab on the stack of your data and the rasterized polygons.

Example, using scot_BNG from example(readOGR) from the rgdal package:

Make some 100x100 rasters containing values 1 to 4:

> r1 = raster(scot_BNG,nrow=100,ncol=100)
> r2 = raster(scot_BNG,nrow=100,ncol=100)
> r1[]=sample(4,ncell(r1),TRUE)
> r2[]=sample(4,ncell(r1),TRUE)

Make a rasterized version of the scottish districts, with the same grid size etc as r1:

> rp = rasterize(scot_BNG, r1)
> plot(rp)

Now make a stack and do a crosstab of layer:

> ct = crosstab(stack(r1,r2,rp))
> head(ct)
  Var1 Var2 Var3 Freq
1    1    1    1    5
2    2    1    1    9
3    3    1    1    7
4    4    1    1    6
5 <NA>    1    1    0

Here Var1 and Var2 are the values from the rasters, and Var3 is the row of the polygon in scot_BNG.

ct has 1425 rows:

> dim(ct)
[1] 1425    4

Which is a row for each of the four values plus NA in the rasters times the number of scottish districts plus 1 for the NA:

> (4+1)*(4+1)*(nrow(scot_BNG)+1)
[1] 1425