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(purrr)

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

uk_bbox_1km <- read_stars(file.path(path_processed, "uk_bbox_1km.tif"))
classes <- 0:7
lcover_files <- file.path(path_cleaned, paste0("landcover_reclass_", classes ,".tif"))
lcover <- map(lcover_files, ~ read_stars(., proxy = TRUE))

Custom functions for computing counts

# counts per axis from original grid to new grid
st_count <- function(src, dest) {
    dd <- dim(dest)
    x_freq <- setNames(rep(0, dd[1]), seq_len(dd[1]))
    y_freq <- setNames(rep(0, dd[2]), seq_len(dd[2]))

    x_seq_g <- st_get_dimension_values(dest, "x", center = FALSE)
    y_seq_g <- st_get_dimension_values(dest, "y", center = FALSE)
    x_seq <- st_get_dimension_values(src, "x", center = NA)
    y_seq <- st_get_dimension_values(src, "y", center = NA)

    x_delta <- st_dimensions(dest)[["x"]]$delta
    y_delta <- st_dimensions(dest)[["y"]]$delta
    x_freq_aux <- table(floor((x_seq - min(x_seq_g)) / x_delta) + 1)
    y_freq_aux <- table(floor((y_seq - max(y_seq_g)) / y_delta) + 1)

    x_freq_aux <- x_freq_aux[names(x_freq_aux) %in% names(x_freq)]
    y_freq_aux <- y_freq_aux[names(y_freq_aux) %in% names(y_freq)]
    x_freq[names(x_freq_aux)] <- x_freq_aux
    y_freq[names(y_freq_aux)] <- y_freq_aux

    return(list(x_freq = x_freq, y_freq = y_freq))
}

# converting average to sum
st_warp_to_sum <- function(src, dest) {
    counts <- st_count(src, dest)
    dest[[1]] <- dest[[1]] * counts[[1]]
    dest[[1]] <- t(t(dest[[1]]) * counts[[2]])
    return(dest)
}

Aggregate land cover at 1 km resolution

lcover_1km <- map(lcover, ~ st_warp(., uk_bbox_1km, use_gdal = TRUE, method = "average"))
lcover_1km <- map2(lcover, lcover_1km, ~ st_warp_to_sum(.x, .y))

Verify the total count per cell is complete

lcover_total <- round(Reduce("+",  lcover_1km))
total_count <-
    prod(sapply(st_dimensions(uk_bbox_1km), function(z) abs(z$delta))) /
    prod(sapply(st_dimensions(lcover[[1]]), function(z) abs(z$delta)))

plot(lcover_total == total_count)

Each cell should sum 1600. In case it is not 1600, we will assign the remaining count to 0 assuming that they do not correspond to any category of interest.

lcover_1km[[1]][[1]] <- lcover_1km[[1]][[1]] + total_count - lcover_total[[1]]
lcover_total <- round(Reduce("+",  lcover_1km))
unique(as.numeric(lcover_total[[1]]))

#> [1] 1600

Save land cover as a raster with bands

lcover_1km <- do.call(c, c(lcover_1km, along = 3))
write_stars(lcover_1km, file.path(path_processed, "uk_1km_landcover.tif"), type = "UInt16")

lcover_1km <- read_stars(file.path(path_processed, "uk_1km_landcover.tif"), proxy = TRUE)
lcover_1km <- st_set_dimensions(lcover_1km, 3, values = classes)

plot_hook = function() plot(uk, border = 'red', add = TRUE)
plot(lcover_1km, hook = plot_hook, main = "Class", breaks = "equal" )

proc.time()

#>    user  system elapsed 
#> 292.542  70.556 371.936