3

I'm trying, in the package terra, to subset a raster by row and column numbers. Apparently that is easy in raster, at least without a geographic extent and crs: Subset a raster using row/column index. But I can't get it to work in terra. There has to be an easy way. terra::subset only selects layers of a raster.

Expecting someone to ask why: I padded a raster with rows and columns before sampling an elevation raster and calculating slope and aspect, which relies on neighboring cells. Now I need to strip off those padded rows and columns.

library(terra)
EXT <- c( -108, -105, 39, 42 )
R <- rast( extent=EXT, ncol=14, nrow=14, crs="epsg:4326" )
R[] <- 1:ncell(R)

# Now try to strip off the outer 2 rows and columns
crop( x=R, y=ext( 3, 12, 3, 12 ) )
# Error: [crop] extents do not overlap

# Normal R-style subsetting also does not work,
# just gives values of that subset
R[ 3:12, 3:12 ]
Jim Worrall
  • 466
  • 3
  • 11

1 Answers1

7

You can subset with drop=FALSE

library(terra)
r <- rast( extent=c( -108, -105, 39, 42 ), ncol=14, nrow=14, crs="epsg:4326" )
values(r) <- 1:ncell(r)

x <- r[ 4:12, 3:10, drop=FALSE]
x
#class       : SpatRaster 
#dimensions  : 9, 8, 1  (nrow, ncol, nlyr)
#resolution  : 0.2142857, 0.2142857  (x, y)
#extent      : -107.5714, -105.8571, 39.42857, 41.35714  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#source      : memory 
#name        : lyr.1 
#min value   :    45 
#max value   :   164

Verify

xy <- xyFromCell(x, cells(x))
range(rowFromY(r, xy[,2]))
#[1]  4 12
range(colFromX(r, xy[,1]))
#[1]  3 10

As to the "why", I would use crop(slope, original) to remove the padded values. That is, the second argument of crop should be your original SpatRaster without the padded cells, or the SpatExtent thereof (ext(orginal))

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63