0

I have two categorical rasters, and I want to get the area of each category on the first further grouped by each category on the second in R, I first tried with zonal stats but you only get the (max, mean, min, count) of the values of the first raster, so now I'm trying to use freq and then multiply by the area the pixel resolution yields, is this ok? Us there a better way to achieve it?

library(raster)
r <- raster(ncols = 5, nrows = 6)
r[] <- sample(rep(1:5, 6), 30)
y <- raster(ncols = 5, nrows = 6)
y[] <- sample(rep(1:3, 10), 30)

df <- data.frame()
for(i in 1:4) { 
  fkwence <- freq(r*(y == i))
  fkwence <- as.data.frame(fkwence)
  fkwence$y_raster_value <- i
  df <- rbind(df, fkwence)
}

and the outcome looks like this:

df
   value count y_raster_value
1      0    20              1
2      1     2              1
3      2     2              1
4      3     2              1
2
  • Is that outcome what you expect?
    – Spacedman
    Dec 28 '19 at 18:56
  • only a part of the outcome, but yes, in the df value would be the r raster's value
    – Elio Diaz
    Dec 30 '19 at 22:10
0

your approach gives correct results but the data frame produced is a bit clumsy to read, I suppose. You would need to subsequently remove the rows with value = 0 because these are artificial values created when using r*(y == i) and that data exists elsewhere in the table (when changing the value of i).

You could use a land use change approach, particularly the change matrix, but using the results to examine the quantity of pixels that exist as value pairs in both rasters. It's a bit more concise as a result, and you can simply sum the freq column by the value of interest, in either raster, for grouping results;

# set up the blank rasters
library(raster)
r <- raster(ncols = 5, nrows = 6)
r[] <- sample(rep(1:5, 6), 30)
y <- raster(ncols = 5, nrows = 6)
y[] <- sample(rep(1:3, 10), 30)

# stack the rasters. Turn into a data frame to get all cellwise value pairs
s <- stack(r,y)
df <- as.data.frame(s)
names(df) <- c("r","y")

# count the rows by unique value pairs
c <- plyr::count(df)

Now you have the quantity of each set of cell values in r per cell value in y

you could quickly visualize in ggplot if you like;

require(ggplot2)
ggplot(data=c,aes(x=r,y=y)) + 
  geom_tile(aes(fill=freq), height=.98, width=.98, na.rm = FALSE) +
  geom_text(aes(x=r,y=y,label=freq),size=6)
2
  • is then just a matter of multiplying by pixel resolution to get area?
    – Elio Diaz
    Jan 3 '20 at 22:54
  • yes that's it, multiply by pixel area (e.g. if they're 500m x 500m, you'll need to multiply each count by 250,000)
    – Sam
    Jan 6 '20 at 7:57

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.