4  Using raster data

This chapter is largely inspired by two presentation; Madelin (2021) and Nowosad (2021); carried out as part of the SIGR2021 thematic school.

4.1 Format of objects SpatRaster

The package terra (Hijmans 2022) allows to handle vector and raster data. To manipulate this spatial data, terra store it in object of type SpatVector and SpatRaster. In this chapter, we focus on the manipulation of raster data (SpatRaster) from functions offered by this package.

An object SpatRaster allows to handle vector and raster data, in one or more layers (variables). This object also stores a number of fundamental parameters that describe it (number of columns, rows, spatial extent, coordinate reference system, etc.).

Source : (Racine 2016)

4.2 Importing and exporting data

The package terra allows importing and exporting raster files. It is based on the GDAL library which makes it possible to read and process a very large number of geographic image formats.

library(terra)

The function rast() allows you to create and/or import raster data. The following lines import the raster file elevation.tif (Tagged Image File Format) into an object of type SpatRaster (default).

elevation <- rast("data_cambodia/elevation.tif") 
elevation
class       : SpatRaster 
dimensions  : 5235, 6458, 1  (nrow, ncol, nlyr)
resolution  : 0.0008333394, 0.0008332568  (x, y)
extent      : 102.2935, 107.6752, 10.33984, 14.70194  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : elevation.tif 
name        : elevation

Modifying the name of the stored variable (altitude).

names(elevation) <- "Altitude"

The function writeRaster() allow you to save an object SpatRaster on your machine, in the format of your choice.

writeRaster(x = elevation, filename = "data_cambodia/new_elevation.tif")

4.3 Displaying a SpatRaster object

The function plot() is use to display an object SpatRaster.

plot(elevation)

A raster always contains numerical data, but it can be both quantitative data and numerically coded qualitative (categorical) data (ex: type of land cover).

Specify the type of data stored with the augment type (type = "continuous" default), to display them correctly.

Import and display of raster containing categorical data: Phnom Penh Land Cover 2019 (land cover types) with a resolution of 1.5 meters:

lulc_2019 <- rast("data_cambodia/lulc_2019.tif")   #Import Phnom Penh landcover 2019, landcover types

The landcover data was produced from SPOT7 satellite image with 1.5 meter spatial resolution. An extraction centered on the municipality of Phnom Penh was then carried out.

plot(lulc_2019, type = "classes")

To display the actual tiles of landcover types you can proceed as follows.

class_name <- c(
  "Roads",
  "Built-up areas",
  "Water Bodies and rivers",
  "Wetlands",
  "Dry bare area",
  "Bare crop fields",
  "Low vegetation areas",
  "High vegetation areas",
  "Forested areas")

class_color <- c("#070401", "#c84639", "#1398eb","#8bc2c2",
                 "#dc7b34", "#a6bd5f","#e8e8e8", "#4fb040", "#35741f")
plot(lulc_2019,
     type = "classes",
     levels = class_name,
     col = class_color,
     plg = list(cex = 0.7),
     mar = c(3.1, 3.1, 2.1, 10)   #The margin are (bottom, left, top, right) respectively
     )

4.4 Change to the study area

4.4.1 (Re)projections

To modify the projection system of a raster, use the function project(). It is then necessary to indicate the method for estimating the new cell values.

Source : Centre Canadien de Télédétection

Four interpolation methods are available:

  • near : nearest neighbor, fast and default method for qualitative data;
  • bilinear : bilinear interpolation. Default method for quantitative data;
  • cubic : cubic interpolation;
  • cubicspline : cubic spline interpolation.
# Re-project data 

elevation_utm = project(x = elevation, y = "EPSG:32648", method = "bilinear")  #from WGS84(EPSG:4326) to UTM zone48N(EPSG:32648) 
lulc_2019_utm = project(x = lulc_2019, y = "EPSG:32648", method = "near") #keep original projection: UTM zone48N(EPSG:32648)

4.4.2 Crop

Clipping a raster to the extent of another object SpatVector or SpatRaster is achievable with the crop().

Source : (Racine 2016)

Import vector data of (municipal divisions) using function vect. This data will be stored in an SpatVector object.

district <- vect("data_cambodia/cambodia.gpkg", layer="district")

Extraction of district boundaries of Thma Bang district (ADM2_PCODE : KH0907).

thma_bang <- subset(district, district$ADM2_PCODE == "KH0907")

Using the function crop(), Both data layers must be in the same projection.

crop_thma_bang <- crop(elevation_utm, thma_bang)

plot(crop_thma_bang)
plot(thma_bang, add=TRUE)

4.4.3 Mask

To display only the values of a raster contained in a polygon, use the function mask().

Source : (Racine 2016)

Creation of a mask on the crop_thma_bang raster to the municipal limits (polygon) of Thma Bang district.

mask_thma_bang <- mask(crop_thma_bang, thma_bang)

plot(mask_thma_bang)
plot(thma_bang, add = TRUE)

4.4.4 Aggregation and disaggregation

Resampling a raster to a different resolution is done in two steps.

1

2

3

Source : (Racine 2016)

Display the resolution of a raster with the function res().

res(elevation_utm)    #check cell size
[1] 91.19475 91.19475

Create a grid with the same extent, then decrease the spatial resolution (larger cells).

elevation_LowerGrid <- elevation_utm
# elevation_HigherGrid  <- elevation_utm

res(elevation_LowerGrid) <- 1000       #cells size = 1000 meter
# res(elevation_HigherGrid) <- 10        #cells size = 10 meter

elevation_LowerGrid
class       : SpatRaster 
dimensions  : 484, 589, 1  (nrow, ncol, nlyr)
resolution  : 1000, 1000  (x, y)
extent      : 203586.3, 792586.3, 1142954, 1626954  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 48N (EPSG:32648)

The function resample() allows to resample the atarting values in the new spatial resolution. Several resampling methods are available (cf. partie 5.4.1).

elevation_LowerGrid <- resample(elevation_utm, 
                                elevation_LowerGrid, 
                                method = "bilinear") 

plot(elevation_LowerGrid, 
     main="Cell size = 1000m\nBilinear resampling method")

4.4.5 Raster fusion

Merge multiple objects SpatRaster into one with merge() or mosaic().

Source : https://desktop.arcgis.com/fr/arcmap/10.3/manage-data/raster-and-images/what-is-a-mosaic.htm

After cutting the elevation raster by the municipal boundary of Thma Bang district (cf partie 5.4.2), we do the same thing for the neighboring municipality of Phnum Kravanh district.

phnum_kravanh <- subset(district, district$ADM2_PCODE == "KH1504")     # Extraction of the municipal boundaries of Phnum Kravanh district

crop_phnum_kravanh <- crop(elevation_utm, phnum_kravanh)             #clipping the elevation raster according to district boundaries

The crop_thma_bang and crop_phnum_kravanh elevation raster overlap spatially:

The difference between the functions merge() and mosiac() relates to values of the overlapping cells. The function mosaic() calculate the average value while merge() holding the value of the object SpaRaster called n the function.

#in this example, merge() and mosaic() give the same result
merge_raster <- merge(crop_thma_bang, crop_phnum_kravanh)   
mosaic_raster <- mosaic(crop_thma_bang, crop_phnum_kravanh)

plot(merge_raster)

# plot(mosaic_raster)

4.4.6 Segregate

Decompose a raster by value (or modality) into different rasterlayers with the function segregate.

lulc_2019_class <- segregate(lulc_2019, keep=TRUE, other=NA)   #creating a raster layer by modality
plot(lulc_2019_class)

4.5 Map Algebra

Map algebra is classified into four groups of operation (Tomlin 1990):

  • Local : operation by cell, on one or more layers;
  • Focal : neighborhood operation (surrounding cells);
  • Zonal : to summarize the matrix values for certain zones, usually irregular;
  • Global : to summarize the matrix values of one or more matrices.

Source : (Li 2009)

4.5.1 Local operations

Source : (Mennis 2015)

4.5.1.1 Value replacement

elevation_utm[elevation_utm[[1]]== -9999] <- NA   #replaces -9999 values with NA

elevation_utm[elevation_utm < 1500]  <- NA        #Replace values < 1500 with NA
elevation_utm[is.na(elevation_utm)] <- 0   #replace NA values with 0

4.5.1.2 Operation on each cell

elevation_1000 <-  elevation_utm + 1000   # Adding 1000 to the value of each cell

elevation_median <-  elevation_utm - global(elevation_utm, median)[[1]]   # Removed median elevation to each cell's value

4.5.1.3 Reclassification

Reclassifying raster values can be used to discretize quantitative data as well as to categorize qualitative categories.

reclassif <- matrix(c(1, 2, 1, 
                      2, 4, 2,
                      4, 6, 3,
                      6, 9, 4), 
                    ncol = 3, byrow = TRUE)

Values between 1 and 2 will be replaced by the value 1.
Values between 3 and 4 will be replaced by the value 2.
Values between 5 and 6 will be replaced by the value 3. Values between 7 and 9 will be replaced by the value 4.

…

reclassif
     [,1] [,2] [,3]
[1,]    1    2    1
[2,]    2    4    2
[3,]    4    6    3
[4,]    6    9    4

The function classify() allows you to perform the reclassification.

lulc_2019_reclass <- classify(lulc_2019, rcl = reclassif)
plot(lulc_2019_reclass, type ="classes")

Display with the official titles and colors of the different categories.

plot(lulc_2019_reclass, 
     type ="classes", 
     levels=c("Urban areas",
              "Water body",
              "Bare areas",
              "Vegetation areas"),
     col=c("#E6004D",
           "#00BFFF",
           "#D3D3D3", 
           "#32CD32"),
     mar=c(3, 1.5, 1, 11))

4.5.1.4 Operation on several layers (ex: NDVI)

It is possible to calculate the value of a cell from its values stored in different layers of an object SpatRaster.

Perhaps the most common example is the calculation of the Normalized Vegetation Index (NDVI). For each cell, a value is calculated from two layers of raster from a multispectral satellite image.

# Import d'une image satellite multispectrale
sentinel2a <- rast("data_cambodia/Sentinel2A.tif")

This multispectral satellite image (10m resolution) dated 25/02/2020, was produced by Sentinel-2 satellite and was retrieved from plateforme Copernicus Open Access Hub. An extraction of Red and near infrared spectral bands, centered on the Phnom Penh city, was then carried out.

plot(sentinel2a)

To lighten the code, we assign the two matrix layers in different SpatRaster objects.

B04_Red <- sentinel2a[[1]]   #spectral band Red

B08_NIR <-sentinel2a[[2]]    #spectral band near infrared

From these two raster objects , we can calculate the normalized vegetation index:

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

raster_NDVI <- (B08_NIR - B04_Red ) / (B08_NIR + B04_Red )

plot(raster_NDVI)

The higher the values (close to 1), the denser the vegetation.

4.5.2 Focal operations

Source : (Mennis 2015)

Focal analysis conisders a cell plus its direct neighbors in contiguous and symmetrical (neighborhood operations). Most often, the value of the output cell is the result of a block of 3 x 3 (odd number) input cells.

The first step is to build a matrix that determines the block of cells that will be considered around each cell.

# 5 x 5 matrix, where each cell has the same weight
mon_focal <- matrix(1, nrow = 5, ncol = 5)
mon_focal
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    1    1    1    1    1
[3,]    1    1    1    1    1
[4,]    1    1    1    1    1
[5,]    1    1    1    1    1

The function focal() Then allows you to perform the desired analysis. For example: calculating the average of the values of all contiguous cells, for each cell in the raster.

elevation_LowerGrid_mean <- focal(elevation_LowerGrid, 
                                  w = mon_focal, 
                                  fun = mean)

4.5.2.1 Focal operations for elevation rasters

The function terrain() allows to perform focal analyzes specific to elevation rasters. Six operations are available:

  • slope = calculation of the slope or degree of inclination of the surface;
  • aspect = calculate slope orientation;
  • roughness = calculate of the variability or irregularity of the elevation;
  • TPI = calculation of the index of topgraphic positions;
  • TRI = elevation variability index calculation;
  • flowdir = calculation of the water flow direction.

Example with calculation of slopes(slope).

#slope calculation
slope <- terrain(elevation_utm, "slope", 
                 neighbors = 8,          #8 (or 4) cells around taken into account
                 unit = "degrees")       #Output unit

plot(slope)                              #Inclination of the slopes, in degrees

4.5.3 Global operations

Source : https://gisgeography.com/map-algebra-global-zonal-focal-local

Global operation are used to summarize the matrix values of one or more matrices.

global(elevation_utm, fun = "mean")  #average values
             mean
Altitude 80.01082
global(elevation_utm, fun = "sd")    #standard deviation
              sd
Altitude 155.885
freq(lulc_2019_reclass)              #frequency
  layer value    count
1     1     1 47485325
2     1     2 13656289
3     1     3 14880961
4     1     4 37194979
table(lulc_2019_reclass[])           #contingency table

       1        2        3        4 
47485325 13656289 14880961 37194979

Statistical representations that summarize matrix information.

hist(elevation_utm)            #histogram
Warning: [hist] a sample of3% of the cells was used

density(elevation_utm)         #density

4.5.4 Zonal operation

Source : (Mennis 2015)

The zonal operation make it possible to summarize the matrix values of certain zones (group of contiguous cells in space or in value).

4.5.4.1 Zonal operation on an extraction

All global operations can be performed on an extraction of cells resulting from the functions crop(), mask(), segregate()…

Example: average elevation for the city of Thma Bang district (cf partie 5.4.3).

# Average value of the "mask" raster over Thma Bang district
global(mask_thma_bang, fun = "mean", na.rm=TRUE)
             mean
Altitude 584.7703

4.5.4.2 Zonal operation from a vector layer

The function extract() allows you to extract and manipulate the values of cells that intersect vector data.

Example from polygons:

# Average elevation for each polygon (district)?
elevation_by_dist <-  extract(elevation_LowerGrid, district, fun=mean)
head(elevation_by_dist, 10)
   ID   Altitude
1   1   8.953352
2   2 196.422240
3   3  23.453937
4   4   3.973118
5   5  29.545801
6   6  41.579593
7   7  50.162749
8   8  85.128777
9   9 269.068091
10 10   8.439041

4.5.4.3 Zonal operation from raster

Zonal operation can be performed by area bounded by the categorical values of a second raster. For this, the two raster must have exaclty the same extent and the same resolution.

#create a second raster with same resolution and extent as "elevation_clip"
elevation_clip <- rast("data_cambodia/elevation_clip.tif")
elevation_clip_utm <- project(x = elevation_clip, y = "EPSG:32648", method = "bilinear")
second_raster_CLC <- rast(elevation_clip_utm)

#resampling of lulc_2019_reclass 
second_raster_CLC <- resample(lulc_2019_reclass, second_raster_CLC, method = "near") 
                               
#added a variable name for the second raster
names(second_raster_CLC) <- "lulc_2019_reclass_resample"

Calculation of the average elevation for the different areas of the second raster.

#average elevation for each area of the "second_raster"
zonal(elevation_clip_utm, second_raster_CLC , "mean", na.rm=TRUE)
  lulc_2019_reclass_resample elevation_clip
1                          1       12.83846
2                          2        8.31809
3                          3       11.41178
4                          4       11.93546

4.6 Transformation and conversion

4.6.1 Rasterization

Convert polygons to raster format.

chamkarmon = subset(district, district$ADM2_PCODE =="KH1201")  
raster_district <- rasterize(x = chamkarmon, y = elevation_clip_utm)
plot(raster_district)

Convert points to raster format

#rasterization of the centroids of the municipalities
raster_dist_centroid <- rasterize(x = centroids(district), 
                                  y = elevation_clip_utm, fun=sum)
plot(raster_dist_centroid, col = "red")
plot(district, add =TRUE)

Convert lines in raster format

#rasterization of municipal boundaries
raster_dist_line <- rasterize(x = as.lines(district), y = elevation_clip_utm, fun=sum)
plot(raster_dist_line)

4.6.2 Vectorisation

Transform a raster to vector polygons.

polygon_elevation <- as.polygons(elevation_clip_utm)
plot(polygon_elevation, y = 1, border="white")

Transform a raster to vector points.

points_elevation <- as.points(elevation_clip_utm)
plot(points_elevation, y = 1, cex = 0.3)

Transform a raster into vector lines.

lines_elevation <- as.lines(elevation_clip_utm)
plot(lines_elevation)

4.6.3 terra, raster, sf, stars…

Reference packages for manipulating spatial data all rely o their own object class. It is sometimes necessary to convert these objects from one class to another class to take advance of all the features offered by these different packages.

Conversion functions for raster data:

FROM/TO raster terra stars
raster rast() st_as_stars()
terra raster() st_as_stars()
stars raster() as(x, ‘Raster’) + rast()

Conversion functions for vector data:

FROM/TO sf sp terra
sf as(x, ‘Spatial’) vect()
sp st_as_sf() vect()
terra st_as_sf() as(x, ‘Spatial’)