Land cover (25 m): Separate classes

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.7.1, GDAL 2.4.0, PROJ 5.2.0

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
#> CRS:            27700

lcover_file <- file.path(path_cleaned, "landcover_reclass.tif")
lcover <- read_stars(lcover_file, proxy = TRUE)

Save binary rasters for each land class

values <- 0:7
lcover_files <- file.path(path_cleaned, paste0("landcover_reclass_", values, ".tif"))
for (i in seq_along(values)) {
    x <- lcover == values[i]
    write_stars(x, lcover_files[i], type = "Byte")
}

#> ====================================================================================================
#> ====================================================================================================
#> ====================================================================================================
#> ====================================================================================================
#> ====================================================================================================
#> ====================================================================================================
#> ====================================================================================================
#> ====================================================================================================

Visualize land cover classes (25 m)

lcovers <- read_stars(lcover_files, proxy = TRUE, along = 3)
lcovers <- st_set_dimensions(lcovers, "new_dim", values = 0:7 )

plot_hook = function() plot(uk, border = 'red', add = TRUE)
plot(lcovers, hook = plot_hook, main = "Class", col = c("white", "black") )

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>     user   system  elapsed 
#> 1005.633  428.649 1435.656