2.7 Raster

The book will focus on regular grids, ie cell with constant size. Other grids exist.

The raster data model consists:

• an header with CRS, extent (cols, rows, cell size), origin (usually lower left)

• a matrix

Raster does not need to store CRS at every cell or at the four corner -> faster processing but a cell can only contains a single value.

2.7.1 R packages for raster data

• {raster} -> {terra}:
  • focus on regular grid,
    
  • one or multi-layered rasters,
    
  • uses C++ or C++ pointers,
    
  • uses its own vector data model,
    
  • rely on built-in function that works on its objects
• {stars}:
  • manage also other grids,
  • raster data cubes,
    
  • stars store value has list of array or file path to large one,
    
  • rely on sf for vector,
    
  • built-in function + methods for existing R functions

In both cases:

• raster should have the same extent

• raster can be read in memory / just the metadata

You can convert from one to an other (st_as_stars() or rast()).

2.7.2 Introduction to terra

Like sp you can still find a lot of raster. terra offer backward compatibility.

terra has the option to divide big raster in smaller chunks.

raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
my_rast = rast(raster_filepath)
class(my_rast)

## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"

my_rast

## class       : SpatRaster 
## dimensions  : 457, 465, 1  (nrow, ncol, nlyr)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : srtm.tif 
## name        : srtm 
## min value   : 1024 
## max value   : 2892

ext(my_rast)

## SpatExtent : -113.239583212784, -112.85208321281, 37.1320834298579, 37.5129167631658 (xmin, xmax, ymin, ymax)

inMemory(my_rast)

## [1] FALSE

# help("terra-package")

2.7.3 Basic map making

plot(my_rast)

# see also terra::plotRGB()
# rasterVis::levelplot()

2.7.4 Raster classes

terra also uses GDAL for reading/writing.

rast() can create raster from scratch:

new_raster = rast(nrows = 6, ncols = 6, 
                  xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
                  vals = 1:36)

SpatRaster can also deals with multiple layers:

multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast = rast(multi_raster_file)
multi_rast

## class       : SpatRaster 
## dimensions  : 1428, 1128, 4  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 12N (EPSG:32612) 
## source      : landsat.tif 
## names       : landsat_1, landsat_2, landsat_3, landsat_4 
## min values  :      7550,      6404,      5678,      5252 
## max values  :     19071,     22051,     25780,     31961

nlyr(multi_rast)

## [1] 4

they can be subset:

multi_rast3 = subset(multi_rast, 3)
multi_rast4 = subset(multi_rast, "landsat_4")

And combine:

multi_rast34 = c(multi_rast3, multi_rast4)

As raster are usualy not stored .rds/.rda are not an option to store them. We can:

• wrap() / unwrap() to save the reference object

• writeRaster()