4.3 Extract information from raster data object

4.3.1 Get CRS

You often need to extract the CRS of a raster object before you interact it with vector data (e.g., extracting values from a raster layer to vector data, or cropping a raster layer to the spatial extent of vector data), which can be done using terra::crs():

terra::crs(IA_cdl_2015_sr)
[1] "BOUNDCRS[\n    SOURCECRS[\n        PROJCRS[\"unnamed\",\n            BASEGEOGCRS[\"GRS 1980(IUGG, 1980)\",\n                DATUM[\"unknown\",\n                    ELLIPSOID[\"GRS80\",6378137,298.257222101,\n                        LENGTHUNIT[\"metre\",1,\n                            ID[\"EPSG\",9001]]]],\n                PRIMEM[\"Greenwich\",0,\n                    ANGLEUNIT[\"degree\",0.0174532925199433,\n                        ID[\"EPSG\",9122]]]],\n            CONVERSION[\"Albers Equal Area\",\n                METHOD[\"Albers Equal Area\",\n                    ID[\"EPSG\",9822]],\n                PARAMETER[\"Latitude of false origin\",23,\n                    ANGLEUNIT[\"degree\",0.0174532925199433],\n                    ID[\"EPSG\",8821]],\n                PARAMETER[\"Longitude of false origin\",-96,\n                    ANGLEUNIT[\"degree\",0.0174532925199433],\n                    ID[\"EPSG\",8822]],\n                PARAMETER[\"Latitude of 1st standard parallel\",29.5,\n                    ANGLEUNIT[\"degree\",0.0174532925199433],\n                    ID[\"EPSG\",8823]],\n                PARAMETER[\"Latitude of 2nd standard parallel\",45.5,\n                    ANGLEUNIT[\"degree\",0.0174532925199433],\n                    ID[\"EPSG\",8824]],\n                PARAMETER[\"Easting at false origin\",0,\n                    LENGTHUNIT[\"metre\",1],\n                    ID[\"EPSG\",8826]],\n                PARAMETER[\"Northing at false origin\",0,\n                    LENGTHUNIT[\"metre\",1],\n                    ID[\"EPSG\",8827]]],\n            CS[Cartesian,2],\n                AXIS[\"easting\",east,\n                    ORDER[1],\n                    LENGTHUNIT[\"metre\",1,\n                        ID[\"EPSG\",9001]]],\n                AXIS[\"northing\",north,\n                    ORDER[2],\n                    LENGTHUNIT[\"metre\",1,\n                        ID[\"EPSG\",9001]]]]],\n    TARGETCRS[\n        GEOGCRS[\"WGS 84\",\n            DATUM[\"World Geodetic System 1984\",\n                ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                    LENGTHUNIT[\"metre\",1]]],\n            PRIMEM[\"Greenwich\",0,\n                ANGLEUNIT[\"degree\",0.0174532925199433]],\n            CS[ellipsoidal,2],\n                AXIS[\"geodetic latitude (Lat)\",north,\n                    ORDER[1],\n                    ANGLEUNIT[\"degree\",0.0174532925199433]],\n                AXIS[\"geodetic longitude (Lon)\",east,\n                    ORDER[2],\n                    ANGLEUNIT[\"degree\",0.0174532925199433]],\n            USAGE[\n                SCOPE[\"Horizontal component of 3D system.\"],\n                AREA[\"World.\"],\n                BBOX[-90,-180,90,180]],\n            ID[\"EPSG\",4326]]],\n    ABRIDGEDTRANSFORMATION[\"Transformation to WGS84\",\n        METHOD[\"Position Vector transformation (geog2D domain)\",\n            ID[\"EPSG\",9606]],\n        PARAMETER[\"X-axis translation\",0,\n            ID[\"EPSG\",8605]],\n        PARAMETER[\"Y-axis translation\",0,\n            ID[\"EPSG\",8606]],\n        PARAMETER[\"Z-axis translation\",0,\n            ID[\"EPSG\",8607]],\n        PARAMETER[\"X-axis rotation\",0,\n            ID[\"EPSG\",8608]],\n        PARAMETER[\"Y-axis rotation\",0,\n            ID[\"EPSG\",8609]],\n        PARAMETER[\"Z-axis rotation\",0,\n            ID[\"EPSG\",8610]],\n        PARAMETER[\"Scale difference\",1,\n            ID[\"EPSG\",8611]]]]"

4.3.2 Subset

You can access specific layers in a multi-layer raster object by indexing:

#--- index ---#
IA_cdl_stack_sr[[2]] # (originally IA_cdl_2016.tif)
class       : SpatRaster 
dimensions  : 11671, 17795, 1  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : -52095, 481755, 1938165, 2288295  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
source      : IA_cdl_2016.tif 
name        : IA_cdl_2016 
min value   :           0 
max value   :         241 

4.3.3 Get cell values

You can access the values stored in a SpatRaster object using the terra::values() function:

#--- terra::values ---#
values_from_rs <- terra::values(IA_cdl_stack_sr)

#--- take a look ---#
head(values_from_rs)
     Layer_1 Layer_1
[1,]       0       0
[2,]       0       0
[3,]       0       0
[4,]       0       0
[5,]       0       0
[6,]       0       0

The returned values come in a matrix form of two columns because we are getting values from a two-layer SpatRaster object (one column for each layer). In general, terra::values() returns a \(X\) by \(n\) matrix, where \(X\) is the number of cells and \(n\) is the number of layers.