List with output from for loop returns empty

I have written a code to obtain crosstab results of a rasterstack for different regions (delimited by a shapefile) covering the raster. However, I am getting an empty list.

This is the function:

transitions <- function(bound, themat) # bound = shapefile # themat = rasterstack
result = vector("list", nrow(bound)) # empty result list
names(result) = bound@data$GEOCODIGO

for (i in 1:nrow(bound)) # this is the number of polygons to iterate through
single <- bound[i,] # selects a single polygon
clip <- mask(crop(themat, single), single) # crops the raster to the polygon boundary

result[i] <- crosstab(clip, digits = 0, long = FALSE, useNA = FALSE)
return(result)




I have tested the steps for the first object in the shapefile/bound outside of the for loop; and it worked well. But I still cannot figure out why I am getting an empty list. Any ideas?

return(result)

for

result[[i]] <- crosstab(...)

[[

1 Answer
1

Example data:

p <- shapefile(system.file("external/lux.shp", package="raster"))
b <- brick(raster(p), nl=2)
values(b) = sample(2, 200, replace=TRUE)


fixed function:

transitions <- function(poly, rast)
result = vector("list", nrow(poly))
for (i in 1:nrow(poly))
clip <- mask(crop(rast, poly[i,]), poly[i,])
result[[i]] <- crosstab(clip, digits = 0, long = FALSE, useNA = FALSE)

return(result)


transitions(p, b)


An alternative would be to use extract

e <- extract(b, p)


To tabulate as in crosstab:

ee <- lapply(e, function(x) aggregate(data.frame(count=rep(1, nrow(x))), data.frame(x), FUN=sum))


To understand that last line, you need to unpack it.

class(e)
#[1] "list"
length(e)
#[1] 12
e[[1]]
# layer.1 layer.2
#[1,] 1 1
#[2,] 1 2
#[3,] 2 2
#[4,] 2 1
#[5,] 2 1
#[6,] 1 2
#[7,] 2 2


e is a list with the same length as the number of polygons (see length(p))

e

length(p)

Let's that the first element and aggregate it to get a table with cases and counts.

x <- e[[1]]

aggregate(data.frame(count=rep(1, nrow(x))), data.frame(x), FUN=sum)
# layer.1 layer.2 count
#1 1 1 1
#2 2 1 2
#3 1 2 2
#4 2 2 2


A similar approach via table (the difference is that you could get Freq values that are zero

as.data.frame(table(x[,1], x[,2]))
# Var1 Var2 Freq
#1 1 1 1
#2 2 1 2
#3 1 2 2
#4 2 2 2


Now wrap the function you like into a lapply

lapply

z <- lapply(e, function(x) aggregate(data.frame(count=rep(1, nrow(x))), data.frame(x), FUN=sum))


And to take it further, bind the data.frames and add an identifier to link the data back to the polygons

y <- do.call(rbind, z,)
y$id <- rep(1:length(z), sapply(z, nrow))

head(y)
# Var1 Var2 Freq id
#1 1 1 1 1
#2 2 1 2 1
#3 1 2 2 1
#4 2 2 2 1
#5 1 1 1 2
#6 2 1 2 2


lapply(e, function(x) aggregate(data.frame(count=rep(1, nrow(x))), data.frame(x), FUN=sum))

Thanks for contributing an answer to Stack Overflow!

But avoid …

To learn more, see our tips on writing great answers.

Required, but never shown

Required, but never shown

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