Search code examples
rgisrasterterra

Subset a raster using row/column index in terra


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 ]

Solution

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