Elevation: aggregate to 1km resolution


Aggregation Copernicus elevation data from 25 m to 1 km resolution. We use the median to aggregate the data to not lose that much information in the borders as happen when using the average method.

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

elev <- read_stars(file.path(path_cleaned, "elev.tif"), proxy = TRUE)
uk_bbox_1km <- read_stars(file.path(path_processed, "uk_bbox_1km.tif"))
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

Visualize current elevation dataset

plot(elev, reset = FALSE, main = "Elevation: 25 m")
plot(st_transform(uk, st_crs(elev)), add = TRUE, border = 2, lwd = 2)

Change resolution in a new CRS

system.time(
    elev_1km <- st_warp(elev, uk_bbox_1km, use_gdal = TRUE, method = "med")
)
#>    user  system elapsed 
#> 125.698  14.539 148.528
write_stars(elev_1km, file.path(path_processed, "uk_1km_elev.tif"))

Visualize new elevation data

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

plot(elev_1km, reset = FALSE, main = "Elevation: 1km")
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 
#> 143.703  16.007 169.894