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