Create UK raster 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

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"), "simple")

#> Reading layer `simple' from data source `/home/rstudio/Documents/Projects/land-suitability/data/cleaned/uk_simple.gpkg' using driver `GPKG'
#> Simple feature collection with 4 features and 10 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -116.1923 ymin: 7054.099 xmax: 655644.8 ymax: 1218625
#> projected CRS:  OSGB 1936 / British National Grid

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

Rasterize UK (1km) using land cover data

Given that we are creating a 1km grid, we make sure that:

• The grid covers all the land cover data by using the bounding box of land cover.
• The reference point (0, 0) of the CRS (27700) is also going to represent a \(1km^2\) cell. We do this by defining the left bottom corner of the bounding box as (-1000, -1000).

# create 1km grid
box <- st_bbox(lcover)
box[c("xmin", "ymin")] <- c(-1000, -1000)
uk_bbox_1km <- st_as_stars(box, dx = 1000, dy = 1000, values = NA_integer_)

# rasterize uk countries
uk_aux <- st_transform(uk, st_crs(lcover))
uk_1km <- st_rasterize(uk_aux, uk_bbox_1km)

# save
write_stars(uk_bbox_1km, file.path(path_processed, "uk_bbox_1km.tif"))
write_stars(uk_1km, file.path(path_processed, "uk_1km_country.tif"))

Visualize 1km rasterized UK

uk_1km <- read_stars(file.path(path_processed, "uk_1km_country.tif"))

plot(uk_1km, reset = FALSE, main = "UK Raster")
plot(st_geometry(uk_aux), add = TRUE, border = 1:4)

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  13.430   0.678  14.055