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)

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"

3 Answers 3


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

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)    
x <- rast(r)
z <- segregate(x, keep=TRUE)


enter image description here

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

#terra 1.6.2
x <- rast(s)
  • Thanks Robert. This gets me a step closer with my own data, but only gives me two values for each plot (like in your reproducible example). However, my values are continuous, so I would need a continuous scalebar.
    – user303287
    Jul 16, 2022 at 9:50
  • 1
    segregate assumes integers. I have expanded the answer to show how you can do this with decimal numbers. Jul 16, 2022 at 16:15
  • I cannot currently install this version: > install.packages("terra", repos = "rspatial.r-universe.dev") Installing package into ‘C:/Users/akuepfer/Documents/R/win-library/4.0’ (as ‘lib’ is unspecified) Package which is only available in source form, and may need compilation of C/C++/Fortran: ‘terra’ These will not be installed > library(terra) terra version 1.4.22
    – user303287
    Jul 17, 2022 at 14:13
  • @user303287: I think that is because your version of R is too old Jul 17, 2022 at 15:10
  • Ah ok. I will update when finished with my deadline :-/, and revisit your solution. Thanks!
    – user303287
    Jul 18, 2022 at 11:00

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.


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]]) )  
  • Thanks for this. I cannot get the first part to work: > sp.terra <- terra::rast(raster::stack(df.sp)) Warning message: NAs introduced by coercion . If I use sp.raster, created from raster(df.sp), then my terra object is without the Unique_id data.
    – user303287
    Jul 14, 2022 at 17:42
  • Alternatively, let's assume we start from a dataframe (which is not sp) with x, y, value, id. How would we go from there using terra? If I understand your reply properly, you are essentially starting from a raster with multiple ids... which is where I need to get to.
    – user303287
    Jul 14, 2022 at 18:12
  • It may be too much to nest both coersions. Try x <- stack(df.sp) then x <- rast(x). Keep in mind that in the raster library raster results in a single band and stack in multiband. For a data.frame you would need to convert to a point feature class then interpolate to a raster. If the points are an array of cells you could use the rasterize function but, need to define an empty raster, defining rows, columns, resolution and bounding coordinates. Jul 14, 2022 at 22:56

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)


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


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


## 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
  • I would double check this approach. Functions such as subset work on data.frame objects but can really muck up sp objects. For sp raster objects the array must remain consistent elsewise offsets occur. This is a questionable approach. If any subset unique ids do not represent the full extent of the raster you will get an error. Jul 14, 2022 at 22:47
  • Thanks for raising possible issues with this approach. It does seem to work for me, and appears to produce exactly the same output for individual unique_ids as it does when I run them one by one. However, just in case, I have amended by answer so that it can also be run on the original dataframe.
    – user303287
    Jul 15, 2022 at 11:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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