example_matrix <- matrix(1:6, nrow = 2, ncol = 3)
example_matrix

##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6

example_matrix[example_matrix >= 3]  <- 3.
example_matrix[example_matrix < 3]  <- 0.
example_matrix

##      [,1] [,2] [,3]
## [1,]    0    3    3
## [2,]    0    3    3

3.1 RasterLayer

A RasterLayer object represents single-layer (variable) raster data. A RasterLayer object always stores a number of fundamental parameters that describe it:

• the number of columns and rows (resolution)
• the spatial extent (cfr. bbox)
• the Coordinate Reference System (CRS)

To create a RasterLayer manually, provide the required minimal information:

# specify the RasterLayer with the following parameters:
# - minimum x coordinate (left border)
# - minimum y coordinate (bottom border)
# - maximum x coordinate (right border)
# - maximum y coordinate (top border)
# - resolution (cell size) in each dimension
r_example <- raster(xmn = 4.42, ymn = 50.9, xmx = 5.4, ymx = 51.4, 
            resolution = c(0.1, 0.1))
r_example

## class       : RasterLayer 
## dimensions  : 5, 10, 50  (nrow, ncol, ncell)
## resolution  : 0.1, 0.1  (x, y)
## extent      : 4.42, 5.42, 50.9, 51.4  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

From the output we learn that with the default parameters, the CRS is defined here in degrees. However, when it would not make sense, the CRS is default put to NA:

raster(ncol = 36, nrow = 18, xmn = -1000, xmx = 1000, ymn = -100, ymx = 900)

## class       : RasterLayer 
## dimensions  : 18, 36, 648  (nrow, ncol, ncell)
## resolution  : 55.55556, 55.55556  (x, y)
## extent      : -1000, 1000, -100, 900  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

Consider the example r_example defined in degrees and extract some additional information:

class(r_example)

## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"

Indeed, we created a single layer raster layer…

res(r_example)

## [1] 0.1 0.1

with a resolution of 0.1 degrees, which we can change:

res(r_example) <- 0.05
r_example

## class       : RasterLayer 
## dimensions  : 10, 20, 200  (nrow, ncol, ncell)
## resolution  : 0.05, 0.05  (x, y)
## extent      : 4.42, 5.42, 50.9, 51.4  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

However, the data currently only defines the skeleton of a raster data set. That is, it knows about its location, resolution, etc., but there are no values associated with it:

hasValues(r_example)

## [1] FALSE

Let’s provide some random numbers from a uniform distribution:

r_example <- setValues(r_example, runif(200, min = 1., max = 10.))
r_example

## class       : RasterLayer 
## dimensions  : 10, 20, 200  (nrow, ncol, ncell)
## resolution  : 0.05, 0.05  (x, y)
## extent      : 4.42, 5.42, 50.9, 51.4  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 1.00307, 9.942678  (min, max)

plot(r_example)

These values can be sliced, just as we would do with vectors and matrices:

values(r_example)[1:10]

##  [1] 8.434002 1.673186 2.903345 8.621475 8.962932 6.403370 4.129534
##  [8] 8.900143 4.838132 6.870831

An important feature of the raster package is the option to work inMemory or not:

inMemory(r_example)

## [1] TRUE

This defines the data itself is stored in the working memory of your computer. For manually created RasterLayers, the data will be in memory. However, when accessing data from files, raster support data access without actually loading all the data in memory.

A useful functionality is the conversion between XYZ columns spatial data representation and a raster representation for points on a regular grid:

xyz_example <- rasterToPoints(r_example)
head(xyz_example)

##          x      y    layer
## [1,] 4.445 51.375 8.434002
## [2,] 4.495 51.375 1.673186
## [3,] 4.545 51.375 2.903345
## [4,] 4.595 51.375 8.621475
## [5,] 4.645 51.375 8.962932
## [6,] 4.695 51.375 6.403370

plot(rasterFromXYZ(xyz_example))

However, in many occasions these XYZ combination are not provided on a regular grid. See the rasterize function for points that are not on a regular grid. An example of points that are not on a regular grid:

longitude <- c(4.4543, 5.02789, 4.94202, 4.49238, 4.49054, 4.54044, 5.95192, 
               4.49496, 4.4958, 5.68327, 4.49054, 4.49054, 4.4958, 4.48938)
latitude <- c(51.33847, 51.24824, 51.24325, 51.26218, 51.25701, 51.27518, 
              51.30803, 51.25803, 51.26106, 51.04185, 51.25701, 51.25701, 
              51.26106, 51.25955)
lonlat <- cbind(longitude, latitude)
head(lonlat)

##      longitude latitude
## [1,]   4.45430 51.33847
## [2,]   5.02789 51.24824
## [3,]   4.94202 51.24325
## [4,]   4.49238 51.26218
## [5,]   4.49054 51.25701
## [6,]   4.54044 51.27518

The output should look like:

3.2 RasterStack and RasterBrick

Single layer raster objects are very common. Still, a collection of rasters with the same spatial extent and resolution is a common case when doing spatial analysis (combining rasters with representing each a variable). In fact, a RasterStack/RasterBrick is a collection of RasterLayer objects with the same spatial extent and resolution.

The main difference between RasterStack and RasterBrick is that a RasterStack is a loose collection of RasterLayer objects that can refer to different files (but must all have the same extent and resolution), whereas a RasterBrick can only point to a single file. A typical example of a RasterBrick object is the representation of a multi-band satellite image.

As an example, create a dummy RasterStack with 3 alternatives of the r_example bundled together:

s_example <- stack(r_example, r_example**2, r_example/2.)
plot(s_example)

When checking the properties

s_example

## class       : RasterStack 
## dimensions  : 10, 20, 200, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 0.05, 0.05  (x, y)
## extent      : 4.42, 5.42, 50.9, 51.4  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## names       :   layer.1,   layer.2,   layer.3 
## min values  : 1.0030703, 1.0061501, 0.5015352 
## max values  :  9.942678, 98.856836,  4.971339

Make sure you notice the additional dimension in the description (nlayers).

3.3 Projection

To transform a RasterLayer to another CRS you can use the function projectRaster. However, whereas the transformation of coordinates (vector formats) is rather trivial, the re-projection of a raster representation can result in a loss of precision and estimates for the values of new cells must be made based on the values in the old cells. If the values are class data, the nearest neighbor is commonly used. Otherwise some sort of interpolation (e.g. bilinear) is required.

An introduction to the CRS class definition to specify the coordinate reference system is provided in the 07-gis-r-vectors.Rmd notebook. The definition of the CRS object is completely similar:

crs_lambert <- CRS("+init=epsg:31370") 

Using the CRS, the re-projection with the projectRaster:

r_reproject <- projectRaster(r_example, crs = crs_lambert)
r_reproject

## class       : RasterLayer 
## dimensions  : 12, 24, 288  (nrow, ncol, ncell)
## resolution  : 3500, 5560  (x, y)
## extent      : 146566.1, 230566.1, 171613.3, 238333.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:31370 +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : 0.9594252, 9.641255  (min, max)

However, using the default parameters, we do not have control on the spatial resolution (dimensions) of the resulting reprojected raster. Setting the resolution provides more control:

r_reproject <- projectRaster(r_example, crs = crs_lambert, 
                             res = 5000)
r_reproject

## class       : RasterLayer 
## dimensions  : 13, 16, 208  (nrow, ncol, ncell)
## resolution  : 5000, 5000  (x, y)
## extent      : 148566.1, 228566.1, 172773.3, 237773.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:31370 +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : 1.346314, 9.209351  (min, max)

But is generally advised to project a raster to another raster object. By providing an existing RasterLayer object, your newly projected data perfectly aligns with it.

# Create a base_layer to project data to
base_raster <- raster(ncol = 20, nrow = 10, 
                      xmn = 153566.1, xmx = 223869.7, 
                      ymn = 176630, ymx = 232773.3)
crs(base_raster) <- crs_lambert
# reproject to the new raster
r_reproject <- projectRaster(r_example, base_raster)
plot(r_reproject)

5.1 Raster algebra

Many generic functions that allow for simple and elegant raster algebra have been implemented for Raster* objects, including the normal algebraic operators such as {}, logical operators such as >, >=, <, ==, ! and functions like abs, round, ceiling, floor, trunc, sqrt, log, log10, exp, cos, sin, atan, tan, max, min, range, prod, sum, any, all. In these functions you can mix raster objects with numbers, as long as the first argument is a raster object.

plot(sqrt(r_example + 10)/2, 
     col = heat.colors(20))

Again, also boolean indexing (conditional replacement) is provided, which supports the reclassification of a raster data set:

r_example[r_example >= 3]  <- 3.
r_example[r_example < 3]  <- 0.
plot(r_example)

The output could look like (depends on your color preferences):

To decide about the color you could use for maps, the colorbrewer website is a great start. You can copy/paste the color codes from the website itself, but the package scales provides the brewer_pal function, providing direct load of the colors as a colormap:

plot(grnt_bodem_reclassified, breaks = c(0, 1, 2, 3), 
     col = scales::brewer_pal(palette = "Greens")(3))

To perform calculations between the individual layers of a RasterBrick or RasterStack, the expression should refer to the individual layers of the object. Referring to individual layers in a RasterBrick or RasterStack object is done by using double square brackets [[]]. As an example, the calculation of the NDVI:

\[NDVI = \frac{NIR - Red}{NIR + Red}\]

ndvi <- (landstat_example[[4]] - landstat_example[[3]]) / (landstat_example[[4]] + landstat_example[[3]])
plot(ndvi)

Although this is a quick way to perform the calculation, directly adding, subtracting, multiplying, etc, the layers of big raster objects is not recommended. When working with big objects, it is advisable to use the calc() function to perform these types of calculations. The reason is that R needs to load all the data first into its internal memory before performing the calculation and then runs everything in one block. It is really easy to run out of memory when doing that. A big advantage of the calc() function is that it has a built-in block processing option for any vectorized function, allowing such calculations to be fully “RAM friendly”. The example below illustrates how to calculate NDVI from the same date set using the calc() function.

## Define the function to calculate NDVI from 
ndvi.calc <- function(x) {
    ndvi <- (x[[4]] - x[[3]]) / (x[[4]] + x[[3]])
    return(ndvi)
}
ndvi2 <- calc(x = landstat_example, fun = ndvi.calc)
plot(ndvi2)

5.2 High-level functions

Several high level functions (i.e. typical GIS software functions) have been implemented for RasterLayer objects. Examples are contour(), focal() (moving window), clump() (detect patches of connected cells), zonal() (zonal statistics), terrain() (slope, aspect and other terrain characteristics from dem), (dis)aggregate (change resolution)… See the help files for more detailed descriptions of each function. (also calc is actually a high level function).

As an example, converting the grnt_bodem raster to a binary mask with 0/1 values can be achieved by the clump function:

plot(clump(grnt_bodem), col = c("darkgreen"))

Another GIS function worthwhile to mention, is the clipping if a raster with a vector layer, which is provided by the function mask:

library(rgdal)
deelbekkens <- readOGR("../data/deelbekkens/Deelbekken.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\stijn_vanhoey\Documents\GitHub\course_gis_scripting\data\deelbekkens\Deelbekken.shp", layer: "Deelbekken"
## with 102 features
## It has 8 fields
## Integer64 fields read as strings:  UIDN OIDN

crs_wgs84 <- CRS("+init=epsg:4326") 
deelbekkens_wgs84 <- spTransform(deelbekkens, crs_wgs84)
r_data_crop <- crop(r_data, deelbekkens_wgs84)
r_data_mask <- mask(r_data_crop, deelbekkens_wgs84)
plot(r_data_mask)