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).


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
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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