Slides : https://nowosad.github.io/SIGR2021/workshop2/workshop2.html#1

srtm_path = system.file("raster/srtm.tif", package = "spDataLarge")
srtm_path
#> [1] "/Users/runner/work/_temp/Library/spDataLarge/raster/srtm.tif"

system.file : nice way to get access to example data contained in packages by getting the path.

R packages are better for storing files than R objects. FIles are mor generic and broader.

Load {terra} and the dataset

library(terra)
#> terra version 1.3.4
srtm = rast(srtm_path)
srtm
#> 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. : +proj=longlat +datum=WGS84 +no_defs 
#> source      : srtm.tif 
#> name        : srtm 
#> min value   : 1024 
#> max value   : 2892

rast() just connect to the file but not loading into memory. It also provide useful information about the dataset.

Looks to alternatives of {terra}: - {raster} : predecessor of {terra}, no new developpement - {stars} : compatible with {sf}, quick, multidimensional, support different types of raster (regular (90% of use cases), rotated, sheared, rectilinear, curvilinear), connect to geospatial servers (works in progress)

Maps

Mapping tools

Maps were detailed earlier, this workshop will use plot() and tmap()

Generally, plot() functions simplifies datasets larger than the image dimensions. It is faster to do that.

With tmap() if there is more than 1 million cells, the data will bo down sampled by tmap()

Map algebra

Four groups - local : computation one cell - focal : cell value and its neigbourg (moving window) - zonal operations summarize raster values for somes zones - global summarize the whole raster (max, min …)

Local operations

raster calculator
srtm2 = srtm * 1000
srtm3 = srtm - globale(srtm, min)[[1]] => 0 -> to max value
srtm3 = srtm - globale(srtm, median)[[1]] => 0 will be the median value, negative value are under the median value
replacing values

srtm_new = srtm
srtm_new[srtm_new < 1500] = NA # values under 1500 will be NA

reclassification
- Reclassification table (2 or 3 columns) : from this value to this value replace by this value
- classify takes reclassification table as argument
operations on many layers (calculating spectral indices, NDVI and such)
- landsat dataset is a simplified dataset from landsat

landsat_path = system.file("raster/landsat.tif", package = "spDataLarge")
landsat = rast(landsat_path)
landsat
#> 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. : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs 
#> source      : landsat.tif 
#> names       : lan_1, lan_2, lan_3, lan_4 
#> min values  :  7550,  6404,  5678,  5252 
#> max values  : 19071, 22051, 25780, 31961

Write a function to compute the ndvi.
use the lapp() function

ndvi_fun = function(nir, red){
  (nir - red) / (nir + red)
}
ndvi = lapp(landsat[[c(4, 3)]], # subset 2 layers (4 and 3) from the dataset
            fun = ndvi_fun) # custom ndvi function

Focal operations

Works on cell neigbourhoood, we need to provide the neighbourghood (window)

The window is provided by a matrix (mandatory type). The function will be applied on the center cell. The matrix can be weighted to have special effect. The matrix needs to be uneven 3x3, 5x5, 7x7. With 4x4 matrix, the function will not be able to determined the center.

Sobel operator : edge detector

The operation starts from the top left to the bottom right, but the order does not change the output.

You can use focal operation on categorical raster : “what’s the most dominant value.”

Zonal operations

also know as zonal statistics. Results are summary table not a raster.

You can provide a function, some are optimised : “mean”, “max”, “min”, “media” note the ".

Global operation

Got a result for the whole raster. global()

Exercices

Find the file path of the nz_elev.tif data using the following code: system.file(“raster/nz_elev.tif”, package = “spDataLarge”). Read this file into R. What is the resolution of the data?

system.file("raster/nz_elev.tif", package = "spDataLarge")
#> [1] "/Users/runner/work/_temp/Library/spDataLarge/raster/nz_elev.tif"

nz_elev <- terra::rast(system.file("raster/nz_elev.tif", package = "spDataLarge"))
nz_elev
#> class       : SpatRaster 
#> dimensions  : 1450, 1115, 1  (nrow, ncol, nlyr)
#> resolution  : 1000, 1000  (x, y)
#> extent      : 995456.5, 2110457, 4741961, 6191961  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
#> source      : nz_elev.tif 
#> name        : elevation 
#> min value   :         0 
#> max value   :  4140.333

The resolution is resolution : 1000, 1000 (x, y)

Create and customize a map based on the nz_elev.tif dataset. Save the result to a .png file.

library(tmap)

nz_map <- tm_shape(nz_elev) +
  #tm_graticules() +
  tm_raster(style = "cont", 
            title = "elevation (m a.s.l)",
            palette = "-Spectral") +
  tm_scale_bar(breaks = c(0, 2, 4),
               text.size = 1) +
  tm_credits("Nicolas Roelandt, 2021") +
  tm_layout(inner.margins = 0,
    main.title = "New Zealand")

tmap_save(nz_map, "vignettes/images/new_zealand_elevation.png")
#> Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown_based_on_GRS80_ellipsoid in Proj4
#> definition
#> stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
#> Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown_based_on_GRS80_ellipsoid in Proj4
#> definition
#> stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
#> Map saved to vignettes/images/new_zealand_elevation.png
#> Resolution: 1841.903 by 2394.263 pixels
#> Size: 6.139675 by 7.980878 inches (300 dpi)
#> Error in dev.off(): QuartzBitmap_Output - unable to open file 'vignettes/images/new_zealand_elevation.png'

Replace values below 0 to 1 (depressions), values between 0 and 300 to 2 (plains), values between 300 and 600 to 3 (hills), and values between 600 and the maximum value to 4 (mountains) in the nz_elev object. Create a new object nz_elev_class.

rcl = matrix(c(-Inf, 0, 1, 
               0, 300, 2, 
               300, 600, 3, 
               600, 4140.333, 4),
             ncol = 3, byrow = TRUE)
rcl
#>      [,1]     [,2] [,3]
#> [1,] -Inf    0.000    1
#> [2,]    0  300.000    2
#> [3,]  300  600.000    3
#> [4,]  600 4140.333    4
nz_elev_class <- classify(nz_elev, rcl = rcl)
nz_elev_class 
#> class       : SpatRaster 
#> dimensions  : 1450, 1115, 1  (nrow, ncol, nlyr)
#> resolution  : 1000, 1000  (x, y)
#> extent      : 995456.5, 2110457, 4741961, 6191961  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
#> source      : memory 
#> name        : elevation 
#> min value   :         2 
#> max value   :         4

Look at the focal function documentation. Apply the Laplacian and the Sobel filters on nz_elev. Compare the results visually.
Calculate an average elevation value in nz_elev for each category in nz_elev_class.
Write a function to calculate the soil-adjusted vegetation index (SAVI; https://en.wikipedia.org/wiki/Soil-adjusted_vegetation_index ). Calculate SAVI for the Landsat 8 images for the area of Zion National Park.

Transformations

Resampling

Reprojecting rasters

different from vector
needs to resample
- create an empty raster template
- to fillthe empty value, use resample()
- parameters : base image, empty raster, function to explain how to resample (method : “bilinear” - good for continues values, “nearest neighbourg” - only one for categorical data)

wgs84 = "EPSG:4326"

nlcd_wgs84 = project(nlcd, wgs84, method = "near")
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'project': object 'nlcd' not found

To create your own projection : projectionwizard.org

Raster-vector interactions

Raster cropping and masking

Cropping : remove data outside the BBox Crop is fast

The mask function replace all values outside the shape with NA. Mask is always used with a cropped image to delete NA data outside the shape

Raster extraction

With points

gps points or anything.

For example, I want elevation of each point using the elevation raster.

terra::extract()

The extract is a dataframe with ID and value.

With polygons

A polygons covers a lotof pixels. Extract each pixel values and give the ID of the polygons.

Useful to see the distribution, do some some statistics.

Rasterization

Transform vector data to raster data

Presence / Absence

rasterize if there is a point in the cell, I get a value of one (TRUE/FALSE)
fun = length : count of point in the cell
field = "vals", fun = sum : sum “vals” values of all points inside the cell

Vectorization

From raster to vector.

Often to create polygons on a categorical raster data

create borders
contours lines (elevation lines, precipitations lines, isolines) : as.contour() : you can provide a number of line to create (nlines) or elevation limits z ,

Conversions

From spatial type/class to another

Output

Save raster file

writeRaster() : depending on the file extension we can have different format With {terra}, geoTIFF are compressed by default.

Raster analysis

Global or local autocorrelation

autocor() : - by default, global autocorrelation - default is ‘Moran’ but you can use ‘Geary’ autocorrelation For elevation, high autocorrelation is expected

For local correlation - provide a window - matrix neighborhood - can be queen (8-n) or rook (4-n) neighborhood - global = FALSE - the output is a raster: find place with high or low autocorrelation

Predictions

{terra} objects can be used on glm, randomForest or prcomp functions.

Interpolations

interpolate()

We have points but we want areas.

Create a model from the points (splines, kriging, IDW)

Geostatistiques with {gstat}

gstat::variogram() distance between values and how values are dissimilear to each other.

Segmentations

{supercells} : super-pixels : automatically find similar areas in visible imagerie

Return an {sf} dataframe with centroids coordinates and mean pixel value of each band of supercell.

Image processing and all things raster

Jakub Nowosad