Slope: compute at 1km resolution

Load packages, read data and source custom scripts

rm(list = ls())
library(stars)

#> Loading required package: abind

#> Loading required package: sf

#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(raster)

#> Loading required package: sp

path_proj <- day2day::git_path()
path_data <- file.path(path_proj, "data")
path_raw <- file.path(path_data, "raw")
path_cleaned <- file.path(path_data, "cleaned")
path_processed <- file.path(path_data, "processed")

uk <- st_read(file.path(path_cleaned, "uk_simple.gpkg"), "union_buffer")

#> Reading layer `union_buffer' from data source `/home/rstudio/Documents/Projects/land-suitability/data/cleaned/uk_simple.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -5115.571 ymin: 2055.298 xmax: 660644.8 ymax: 1223625
#> projected CRS:  OSGB 1936 / British National Grid

elev_1km <- read_stars(file.path(path_processed, "uk_1km_elev.tif"), proxy = FALSE)
elev_crs <- st_crs(elev_1km)

Fill missing values

Filling missing values (bathymetry) with zero values to not lose information when computing the slope.

id_na <- is.na(elev_1km[[1]])
elev_1km[[1]][id_na] <- 0

plot(elev_1km, reset = FALSE, main = "Auxiliary elevation")
plot(st_transform(uk, st_crs(elev_1km)), add = TRUE, border = 2, lwd = 2)

Compute slope with 4 and 8 neighbors

# computing slope with raster package
elev_1km <- as(elev_1km, "Raster")

#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS definition,
#>  but +towgs84= values preserved

slope_nb4 <- terrain(elev_1km, "slope", "degrees", neighbors = 4)
slope_nb8 <- terrain(elev_1km, "slope", "degrees", neighbors = 8)

# convert to stars
slope_nb4 <- st_set_crs(st_as_stars(slope_nb4), elev_crs)
slope_nb8 <- st_set_crs(st_as_stars(slope_nb8), elev_crs)

# remove artificial bathymetry slope
slope_nb4[[1]][id_na] <- NA
slope_nb8[[1]][id_na] <- NA

write_stars(slope_nb4, file.path(path_processed, "uk_1km_slope_nb4.tif"))
write_stars(slope_nb8, file.path(path_processed, "uk_1km_slope_nb8.tif"))

Visualize slope data

slope_nb4 <- read_stars(file.path(path_processed, "uk_1km_slope_nb4.tif"), proxy = TRUE)
slope_nb8 <- read_stars(file.path(path_processed, "uk_1km_slope_nb8.tif"), proxy = TRUE)

plot(slope_nb4, reset = FALSE, main = "Slope: 4 neighbors", col = hcl.colors(10))
plot(st_transform(uk, st_crs(elev_1km)), add = TRUE, border = 2, lwd = 2)

plot(slope_nb8, reset = FALSE, main = "Slope: 8 neighbors", col = hcl.colors(10))
plot(st_transform(uk, st_crs(elev_1km)), add = TRUE, border = 2, lwd = 2)

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  24.215   0.873  25.226