4.1 Raster data object classes

4.1.1 `raster` package: `RasterLayer`, `RasterStack`, and `RasterBrick`

Let’s start with taking a look at raster data. We will use the CDL data for Iowa in 2015.

class      : RasterLayer 
dimensions : 11671, 17795, 207685445  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : -52095, 481755, 1938165, 2288295  (xmin, xmax, ymin, ymax)
crs        : +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_2015.tif 
names      : IA_cdl_2015 
values     : 0, 229  (min, max)

Evaluating the imported raster object provides you with information about the raster data, such as dimensions (number of cells, number of columns, number of cells), spatial resolution (30 meter by 30 meter for this raster data), extent, CRS and the minimum and maximum values recorded in this raster layer. The class of the downloaded data is RasterLayer, which is a raster data class defined by the raster package. A RasterLayer consists of only one layer, meaning that only a single variable is associated with the cells (here it is land use category code in integer).

You can stack multiple raster layers of the same spatial resolution and extent to create a RasterStack using raster::stack() or RasterBrick using raster::brick(). Often times, processing a multi-layer object has computational advantages over processing multiple single-layer one by one⁶⁷.

To create a RasterStack and RasterBrick, let’s load the CDL data for IA in 2016 and stack it with the 2015 data.

class      : RasterStack 
dimensions : 11671, 17795, 207685445, 2  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : -52095, 481755, 1938165, 2288295  (xmin, xmax, ymin, ymax)
crs        : +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 
names      : IA_cdl_2015, IA_cdl_2016 
min values :           0,           0 
max values :         229,         241

IA_cdl_stack is of class RasterStack, and it has two layers of variables: CDL for 2015 and 2016. You can make it a RasterBrick using raster::brick():

#--- stack the two ---#
IA_cdl_brick <- brick(IA_cdl_stack)

#--- or this works as well ---#
# IA_cdl_brick <- brick(IA_cdl_2015, IA_cdl_2016)

#--- take a look ---#
IA_cdl_brick

class      : RasterBrick 
dimensions : 11671, 17795, 207685445, 2  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : -52095, 481755, 1938165, 2288295  (xmin, xmax, ymin, ymax)
crs        : +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     : r_tmp_2022-04-12_214849_67339_89224.grd 
names      : IA_cdl_2015, IA_cdl_2016 
min values :           0,           0 
max values :         229,         241

You probably noticed that it took some time to create the RasterBrick object⁶⁸. While spatial operations on RasterBrick are supposedly faster than RasterStack, the time to create a RasterBrick object itself is often long enough to kill the speed advantage entirely⁶⁹. Often, the three raster object types are collectively referred to as Raster\(^*\) objects for shorthand in the documentation of the raster and other related packages.

4.1.2 `terra` package: `SpatRaster`

terra package has only one object class for raster data, SpatRaster and no distinctions between one-layer and multi-layer rasters is necessary. Let’s first convert a RasterLayer to a SpatRaster using terra::rast() function.

#--- convert to a SpatRaster ---#
IA_cdl_2015_sr <- rast(IA_cdl_2015)

#--- take a look ---#
IA_cdl_2015_sr

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_2015.tif 
name        : IA_cdl_2015 
min value   :           0 
max value   :         229

You can see that the number of layers (nlyr in dimensions) is \(1\) because the original object is a RasterLayer, which by definition has only one layer. Now, let’s convert a RasterStack to a SpatRaster using terra::rast().

#--- convert to a SpatRaster ---#
IA_cdl_stack_sr <- rast(IA_cdl_stack)

#--- take a look ---#
IA_cdl_stack_sr

class       : SpatRaster 
dimensions  : 11671, 17795, 2  (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 
sources     : IA_cdl_2015.tif  
              IA_cdl_2016.tif  
names       : IA_cdl_2015, IA_cdl_2016 
min values  :           0,           0 
max values  :         229,         241

Again, it is a SpatRaster, and you now see that the number of layers is 2. We just confirmed that terra has only one class for raster data whether it is single-layer or multiple-layer ones.

In order to make multi-layer SpatRaster from multiple single-layer SpatRaster you can just use c() like below:

#--- create a single-layer SpatRaster ---#
IA_cdl_2016_sr <- rast(IA_cdl_2016)

#--- concatenate ---#
(
  IA_cdl_ml_sr <- c(IA_cdl_2015_sr, IA_cdl_2016_sr)
)

(
  IA_cdl_ml_sr <- c(IA_cdl_2015_sr, IA_cdl_2016_sr)
)

class       : SpatRaster 
dimensions  : 11671, 17795, 2  (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 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
sources     : IA_cdl_2015.tif  
              IA_cdl_2016.tif  
names       : Layer_1, Layer_1 
min values  :       0,       0 
max values  :     229,     241

4.1.3 Converting a `SpatRaster` object to a `Raster`\(^*\) object.

You can convert a SpatRaster object to a Raster\(^*\) object using raster(), stack(), and brick(). Keep in mind that if you use rater() even though SpatRaster has multiple layers, the resulting RasterLayer object has only the first of the multiple layers.

#--- RasterLayer (only 1st layer) ---#
IA_cdl_stack_sr %>% raster()

class      : RasterLayer 
dimensions : 11671, 17795, 207685445  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : -52095, 481755, 1938165, 2288295  (xmin, xmax, ymin, ymax)
crs        : +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_2015.tif 
names      : IA_cdl_2015 
values     : 0, 229  (min, max)

#--- RasterLayer ---#
IA_cdl_stack_sr %>% stack()

class      : RasterStack 
dimensions : 11671, 17795, 207685445, 2  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : -52095, 481755, 1938165, 2288295  (xmin, xmax, ymin, ymax)
crs        : +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 
names      : IA_cdl_2015, IA_cdl_2016 
min values :           0,           0 
max values :         229,         241

#--- RasterLayer (this takes some time) ---#
IA_cdl_stack_sr %>% brick()

class      : RasterStack 
dimensions : 11671, 17795, 207685445, 2  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : -52095, 481755, 1938165, 2288295  (xmin, xmax, ymin, ymax)
crs        : +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 
names      : IA_cdl_2015, IA_cdl_2016 
min values :           0,           0 
max values :         229,         241

Or, you can just use as(SpatRast, "Raster") like below:

as(IA_cdl_stack_sr, "Raster")

class      : RasterStack 
dimensions : 11671, 17795, 207685445, 2  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : -52095, 481755, 1938165, 2288295  (xmin, xmax, ymin, ymax)
crs        : +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 
names      : IA_cdl_2015, IA_cdl_2016 
min values :           0,           0 
max values :         229,         241

This works for any Raster\(^*\) object and you do not have to pick the right function like above.

4.1.4 Vector data in the `terra` package

terra package has its own class for vector data, called SpatVector. While we do not use any of the vector data functionality provided by the terra package, we learn how to convert an sf object to SpatVector because terra functions do not support sf as of now (this will likely be resolved very soon). We will see some use cases of this conversion in the next chapter when we learn raster value extractions for vector data using terra::extract().

As an example, let’s use Illinois county border data.

Simple feature collection with 102 features and 17 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -91.51308 ymin: 36.9703 xmax: -87.01993 ymax: 42.50848
Geodetic CRS:  NAD83
First 10 features:
    STATEFP COUNTYFP COUNTYNS GEOID    NAME       NAMELSAD LSAD CLASSFP MTFCC
86       17      067 00424235 17067 Hancock Hancock County   06      H1 G4020
93       17      025 00424214 17025    Clay    Clay County   06      H1 G4020
132      17      185 00424293 17185  Wabash  Wabash County   06      H1 G4020
149      17      113 01784833 17113  McLean  McLean County   06      H1 G4020
159      17      005 00424204 17005    Bond    Bond County   06      H1 G4020
160      17      009 00424206 17009   Brown   Brown County   06      H1 G4020
214      17      083 00424243 17083  Jersey  Jersey County   06      H1 G4020
255      17      147 00424275 17147   Piatt   Piatt County   06      H1 G4020
267      17      151 00424277 17151    Pope    Pope County   06      H1 G4020
304      17      011 00424207 17011  Bureau  Bureau County   06      H1 G4020
    CSAFP CBSAFP METDIVFP FUNCSTAT      ALAND   AWATER    INTPTLAT     INTPTLON
86    161  22800     <NA>        A 2055798688 53563362 +40.4013180 -091.1688008
93   <NA>   <NA>     <NA>        A 1213022841  3064720 +38.7468187 -088.4823254
132  <NA>   <NA>     <NA>        A  578403995 10973569 +38.4458209 -087.8391674
149   145  14010     <NA>        A 3064597287  7804856 +40.4945594 -088.8445391
159   476  41180     <NA>        A  985065096  6462629 +38.8859240 -089.4365916
160  <NA>   <NA>     <NA>        A  791828626  4144346 +39.9620694 -090.7503095
214   476  41180     <NA>        A  957415216 20333974 +39.0801945 -090.3613850
255  <NA>  16580     <NA>        A 1137492089   754122 +40.0090327 -088.5923546
267  <NA>   <NA>     <NA>        A  955362037 14662240 +37.4171687 -088.5423737
304   176  36837     <NA>        A 2250884551 11523914 +41.4013043 -089.5283772
                          geometry
86  MULTIPOLYGON (((-91.37421 4...
93  MULTIPOLYGON (((-88.69517 3...
132 MULTIPOLYGON (((-87.9446 38...
149 MULTIPOLYGON (((-89.2665 40...
159 MULTIPOLYGON (((-89.36179 3...
160 MULTIPOLYGON (((-90.91469 4...
214 MULTIPOLYGON (((-90.59216 3...
255 MULTIPOLYGON (((-88.74516 4...
267 MULTIPOLYGON (((-88.70963 3...
304 MULTIPOLYGON (((-89.85691 4...

You can convert an sf object to SpatVector object using terra::vect().

(
  IL_county_sv <- vect(IL_county)
)

 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 102, 17  (geometries, attributes)
 extent      : -91.51308, -87.01993, 36.9703, 42.50848  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat NAD83 (EPSG:4269) 
 names       : STATEFP COUNTYFP COUNTYNS GEOID    NAME       NAMELSAD  LSAD
 type        :   <chr>    <chr>    <chr> <chr>   <chr>          <chr> <chr>
 values      :      17      067 00424235 17067 Hancock Hancock County    06
                    17      025 00424214 17025    Clay    Clay County    06
                    17      185 00424293 17185  Wabash  Wabash County    06
 CLASSFP MTFCC CSAFP (and 7 more)
   <chr> <chr> <chr>             
      H1 G4020   161             
      H1 G4020    NA             
      H1 G4020    NA

You will see this in Chapter 5 where we learn how to extract values from a raster layer for a vector data.↩︎
Read here for the difference between RasterStack and RasterBrick↩︎
We will see this in Chapter , where we compare the speed of data extraction from RasterStack and RasterBrick objects.↩︎