Create multiple rasters from single spatialpixeldataframe, by Unique id

Question

I have a SpatialPixelDataFrame (spdf) which contains gridded values (UD) for different Unique id's (n = 21).

I would like to create an individual raster for each Unique_id, with value (UD) as the cell value.

I can create a single raster using the following:

single.ras <- raster(spdf)

Is there a way to group the data by Unique_id and produce multiple raster files (e.g. in the form of a Raster stack?)

For example, for an estUDm object, I have the following code which creates a stack of rasters, one for each Unique id. However, the difference here is that the individual Unique_id "objects" are already inherent in the estUDm object.

multi.ras <- stack(lapply(estUDm.data, raster))

I have tried to adapt a code that uses a loop, but this just produces the same raster 16 (? why 16?) times. So something is not correct (I know next to nothing about loops).

## Loop
datagroup <- unique(df.sp$Unique_id)

result <- list()
for (i in 1:length(datagroup)) {
  result[[i]] <- raster(df.sp)
}  
s <- stack(result)
plot(s)

The data structure of my spdf is as follows; apologies I do not know how to recreate such data.

> str(spdf)
Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
  ..@ data       : tibble [623 x 2] (S3: tbl_df/tbl/data.frame)
  .. ..$ UD       : num [1:623] 0.1064 0.0997 0.0937 0.0346 0.0414 ...
  .. ..$ Unique_id: chr [1:623] "Fixed gear_2013_B" "Fixed gear_2013_B" "Fixed gear_2013_B" "Fixed gear_2013_B" ...
  ..@ coords.nrs : num(0) 
  ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
  .. .. ..@ cellcentre.offset: Named num [1:2] -69.2 -55.8
  .. .. .. ..- attr(*, "names")= chr [1:2] "Lon" "Lat"
  .. .. ..@ cellsize         : Named num [1:2] 0.5 0.5
  .. .. .. ..- attr(*, "names")= chr [1:2] "Lon" "Lat"
  .. .. ..@ cells.dim        : Named int [1:2] 27 14
  .. .. .. ..- attr(*, "names")= chr [1:2] "Lon" "Lat"
  ..@ grid.index : int [1:623] 247 220 166 139 112 85 248 248 221 194 ...
  ..@ coords     : num [1:623, 1:2] -67.8 -67.8 -67.8 -67.8 -67.8 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:623] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "Lon" "Lat"
  ..@ bbox       : num [1:2, 1:2] -69.5 -56 -56 -49
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "Lon" "Lat"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"

Answer 1

You can use terra::segregate to create a SpatRaster with a layer for each class (unique ID) that is present in a single layer input SpatRaster.

Example data

library(raster)
r <- raster(system.file("external/test.grd", package="raster")) |> cut(5)
s <- as(r, "SpatialPixelsDataFrame")

Create a RasterLayer and from that a SpatRaster. Then use terra::segregate. The default value assigned to the other cells that are not NA is zero, but you can change that with argument other.

r <- raster(s)    
library(terra)
x <- rast(r)
z <- segregate(x, keep=TRUE)

plot(z)

There is also a short algebraic approach that gives you a Boolean layer for each ID:

 y <- x == 1:5

If your values are decimal numbers, not integers, for example

 xx <- x / 3

you can do

 u <- unlist(unique(xx))
 yy <- (xx == u) * u

Since "terra" version 1.6.2 you can also use segregate with decimal numbers

zz <- segregate(xx, keep=TRUE)
zz
#class       : SpatRaster 
#dimensions  : 115, 80, 5  (nrow, ncol, nlyr)
#resolution  : 40, 40  (x, y)
#extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#names       :  0.333333,  0.666667,         1,  1.333333,  1.666667 
#min values  :         0,         0,         0,         0,         0 
#max values  : 0.3333333, 0.6666667, 1.0000000, 1.3333333, 1.6666667 

Jeffrey Evens noted that "there are no direct coercion methods from sp to terra". That has now been fixed for all "sp" classes. With "terra" version 1.6.2 you can skip the "raster" bit, and get a SpatRaster from a SpatialPixels object like this

library(terra)
#terra 1.6.2
x <- rast(s)

Answer 2

I would first recommend getting out of an sp class and second, use terra in leu of raster. This moves you into modern libraries (sp and raster are on a depreciation schedule with sf and terra as replacements) and opens up the use of terra::ifel, which simplifies things considerably.

There are no direct coercion methods from sp to terra so, you can use the raster library as a stepping stone, coercing sp to a raster stack then terra rast using terra::rast(raster::stack(sp.raster)).

Here is a worked example of creating separate rasters based on an ID raster.

Add library and create a random raster (r) and an id raster ids. We then combine them into a single raster stack. THis is not necessary but, emulates what you would get converting an sp multiband object to a raster stack to a terra rast object.

library(terra)

r <- rast(nrow=20, ncol=20)
  r[] <- runif(ncell(r))
ids <- r
  ids[] <- sample(1:2, ncell(r), replace=TRUE)
( r <- c(ids, r) )

Now, we iterate through the unique ID's, in the ids raster, to subset into separate rasters. The double bracket denote which raster in the stack we are operating on. The call to unique(r[[1]])[,1] returns the unique values of the raster but, since it returns a data.frame we have to subset the first column thus, the [,1] index. We then use ifel (as special case of ifelse for use on terra rasters). We are saying if ids == i, return associated cells in r else return NA. rids <- lapply(unique(r[[1]])[,1], function(i) { ifel(r[[1]] == i, r[[2]], NA) })
plot( rids <- c(rids[[1]], rids[[2]]) )

This can easily be expanded to include multiple raster subsets.

r <- c(r, r[[2]]*10) # add 3rd raster

# Note; r[[2:3]] subsetting 2nd and 3rd rasters 
rids <- lapply(unique(r[[1]])[,1], function(i) { 
  ifel(r[[1]] == i, r[[2:3]], NA)
})  
plot( rids <- c(rids[[1]], rids[[2]]) )  

Answer 3

The terra route is probably the best way to go given the changes in packages. However, if you need a fix with the raster package, a friend has fixed my half baked loop, so this works using sp and raster package:

datagroup <- unique(spdf$Unique_id) ## where spdf is a specialpixeldataframe
result <- list()
for (i in 1:length(datagroup)) {
  df.unique <- subset(spdf, Unique_id == datagroup[i])
  result[[i]] <- raster(df.unique) ## this assumes that the values for the cells are in your first column. if the values are in e.g. column 4, write raster(df.unique[4])
}  

names(result) <- datagroup ## to make your life easier - this will identify each of your raster files by its Unique_id.
s <- stack(result)

plot(s)

writeRaster(stack(s), names(s), bylayer=TRUE, format='GTiff')#stackOpen(stackfile)

EDIT:

if subsetting of an sp object causes issues, the following code could be used on a dataframe.

library(rasterize)
library(raster)

## produce an empty raster
r = raster(ext=extent(c(-69.5, -56, -56, -49)), res=c(0.5, 0.5)) 
r[] <- 1
projection(r) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
crs(r) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

## Convert dataframe to multiple rasters
datagroup <- unique(df$Unique_id) ## choose your id
result <- list()
for (i in 1:length(datagroup)) {
  df.unique <- subset(df, Unique_id == datagroup[i]) ## subset by unique_id
  result[[i]] <- rasterize(x = df.unique[, 2:3], y= r, field =df.unique[,5], fun = sum ) ## x = lon,lat; y = empty raster object; field = values
result[[i]][is.na(result[[i]][])] <- 0  ## convert NAs to 0
  
}