Grid

Summary

This tutorial explore how to make a grid and transform observations into a grid with terra package:

import a multi-layer rasters with terra::rast()
extract information from rasters on points or polygons with terra::extract()

The ecologist mind

If we want to make some modeling, it can be needed to regroup observations into a grid.

Setup

Follow the setup instructions if you haven’t followed previous tutorials

If haven’t done it already, please follow the setup instructions.

Let’s start with loading the required packages.

suppressPackageStartupMessages({
  library(mapview)
  library(here)
  library(terra)
  # library(sf)
  # library(exactextractr)
})

Now load all the datasets attached to the observations of otters recorded in 2021 within a 50km buffer from Montpellier, France.

pt_otter <- vect(here("data", "gbif_otter_2021_mpl50km.gpkg"))
river <- vect(here("data", "BDCARTO-River_mpl50km.gpkg"))
landuse <- vect(here("data", "BDCARTO-LULC_mpl50km.shp"))
bdalti <- rast(here("data", "BDALTI_mpl50km.tif"))
temperature <- rast(here("data", "CHELSA_monthly_tas_2015_2021.tif"))

Create a grid

Key issue : which resolution? which projection? resolution depends on your dataset, extent, and ecological question projection depends on the dataset

Let’s start from our observation.

pt_2154 <- project(pt_otter, "EPSG:2154")
res <- 5000 # 1km grid
grid <- rast(pt_2154, res = res)
grid

class       : SpatRaster 
size        : 17, 19, 1  (nrow, ncol, nlyr)
resolution  : 5000, 5000  (x, y)
extent      : 723524.6, 818524.6, 6245001, 6330001  (xmin, xmax, ymin, ymax)
coord. ref. : RGF93 v1 / Lambert-93 (EPSG:2154)

Points to grid

the argument fun=count set that it calculate the number of points the argument background=0 set that f there is no point, the vale will be 0, instead of NA by default

grid_nobs <- rasterize(pt_2154, grid, fun = "count", background = 0)

plot(grid_nobs)

Lines to grid

length of rivers

river_2154 <- project(river, "EPSG:2154")

grid_river <- rasterizeGeom(river_2154, grid, fun = "length")
plot(grid_river)

Polygon to grid

percentage of forest

# select only the forest
forest <- landuse[landuse$nature == "Forêt"]
# project the polygons
forest_2154 <- project(forest, "EPSG:2154")
# rasterize
grid_forest <- rasterize(forest_2154, grid, cover = TRUE, background = 0)

plot(grid_forest)

Resample raster

Elevation

grid_alti <- resample(bdalti, grid, "bilinear")
plot(grid_alti)

Temperature

#calculate the average temperature in 2021
in_2021 <- grep("2021", names(temperature))
monthly_temp <- subset(temperature, in_2021)
annual_temp <- mean(monthly_temp, na.rm = TRUE)
# project (that's one of the only occasion when projection of raster)
temp_2154 <- project(annual_temp, "EPSG:2154")
grid_temp <- resample(temp_2154, grid, "bilinear")
plot(grid_temp)

Merge all together

Ignoring all concept of spatial bias or anything else

out <- c(grid_nobs, grid_river, grid_forest, grid_alti, grid_temp)

names(out) <- c("nobs", "river_m", "forest", "elevation_m", "temperature_K")

plot(out)

france_border <- geodata::gadm("FRA", level = 0, path = here("data"))
# project the coast in Lambert-93
france_2154 <- project(france_border, crs(bdalti))
# mask (=set to NA) the pixels that are not in the polygon
out_masked <- mask(out, france_2154)
plot(out_masked)

df <- data.frame(out_masked)

m1 <- lm(nobs ~ river_m + forest + elevation_m + temperature_K, data = df)
anova(m1)

Analysis of Variance Table

Response: nobs
               Df  Sum Sq Mean Sq F value   Pr(>F)   
river_m         1   0.360  0.3602  0.5795 0.447163   
forest          1   0.068  0.0682  0.1097 0.740722   
elevation_m     1   0.623  0.6229  1.0021 0.317683   
temperature_K   1   5.463  5.4635  8.7896 0.003294 **
Residuals     276 171.556  0.6216                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Your turn

Does the spatial resolution has en effect on the results? Repeat the same process with a grid of 500m, 1km, or 10km resolution and compare the results.