0

I have a raster file of land cover which I have reduced to contain just tree cover cells. I have used clump in the raster package to clump() together contiguous areas of forest. This gives all cells touching one another the same ID because they are part of the same patch.
I then want to figure out the PatchStat() for each clump, which I do by converting my clump raster to as.matrix. I tried to get the PatchStat() to do this to the raster but it would only work if it was in a matrix.

I now want to make a raster with the patch stat output, namely "perim.area.ratio". So each cell that corresponds to clump 1, will get the perim.area.ratio value that corresponds to clump 1. To do this I made a data.frame() from my clump raster that had: lon, lat, layer(clumpID), cellID.
I tried to merge my clump raster data.frame with the PatchStat output using layer and patchID. However, an error occurs:

Error in fix.by(by.x, x) : 'by' must specify valid column(s).

Any ideas how I could do this another way, or why these columns aren't valid? Code below.

clump <- raster(file.choose())
library(SDMTools)
clumpval <- rasterToPoints(clump)
clumpcell <- cellFromXY(clump, clumpval[, c('x', 'y')] )
clumpdf <- data.frame(clumpval, clumpcell)
ps.data <- PatchStat(as.matrix(clump))
merged.data.all <- merge(clumpdf, ps.data1, by=c("layer", "patchID"))
jbaums
  • 27,115
  • 5
  • 79
  • 119
Adam
  • 1,147
  • 3
  • 15
  • 23
  • The help file (http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=SDMTools:PatchStat) says that: "a matrix of data with individual patches identified as with ConnCompLabel; The matrix can be a raster of class 'asc' (adehabitat package), 'RasterLayer' (raster package) or 'SpatialGridDataFrame' (sp package)". If you are working with RasterLayer (if you're using `raster` you probably are), you could avoid having to convert to matrix and back. – Roman Luštrik Mar 15 '12 at 08:19
  • Hi @RomanLustirk, thanks for the reply. Yes I am working with a 'RasterLayer', when I check the properties of my raster it says it at the top. However, it does't work if I just put in the raster, I get an error message... Not sure why but couldn't figure that out either! Cheers, Adam – Adam Mar 15 '12 at 09:59
  • What is the error message? `r <- raster(ncols=12, nrows=12); r[] <- round(runif(ncell(r))*0.6 ); rc <- clump(r); PatchStat(rc)` works for me. – jbaums Mar 15 '12 at 13:01
  • G'day @jbaums this is the error i get: Error in if (m[ii] == 0) minp[ii] = 4 * n[ii] : missing value where TRUE/FALSE needed – Adam Mar 20 '12 at 02:19
  • I can't reproduce your 'missing value where TRUE/FALSE needed' error, but have updated my solution below to explain the reason for your primary problem. – jbaums Mar 21 '12 at 05:57

1 Answers1

2

The way you've coded it, the merge function is expecting both data frames to have both a 'layer' field and a 'patchID' column, whereas in fact you intend to map the layer column of clumpdf to the patchID column of ps.data. You need to use by.x and by.y arguments.

The correct call would be: 

merged.data.all <- merge(clumpdf, ps.data, by.x='layer', by.y="patchID")

However, there's another simple way to assign cells their clump's perim.area.ratio:

library(raster)
library(SDMTools)

# create a random raster
r <- raster(ncols=200, nrows=200)
r[] <- rbinom(ncell(r), 1, 0.5)

# clump it
rc <- clump(r)

# get patch stats
p <- PatchStat(rc)

# Replace each non-NA value of rc with the corresponding clump perim.area.ratio.
not.na <- Which(!is.na(rc), cells=TRUE)
rc[not.na] <- sapply(rc[not.na], function(x) {
  p[p$patchID==x, 'perim.area.ratio']
})

To break it down for you (in case you're not particularly familiar with apply functions), this last bit firstly identifies the cell indices of all cells with non-NA values, and assigns this vector to the object not.na. The sapply function then assigns each value of not.na in turn to x, and performs the stuff between the curly braces (which in this case just returns the perim.area.ratio value found in the row of p that has patchID equal to x). The sapply function returns a vector of these perim.area.ratio values, which is then assigned to the non-NA cells of rc. In essence, it's a find & replace operation, where patch numbers are replaced by their corresponding perim.area.ratios.

I should mention that this may not work if you have a very large grid.

jbaums
  • 27,115
  • 5
  • 79
  • 119