Elevation: merge and crop elevation data for UK

Load packages, read data and source custom scripts

rm(list = ls())
library(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag

#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

library(ggplot2)
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_dem <- file.path(path_raw, "elevation")

tile3030 <- read_stars(file.path(path_dem, "eu_dem_v11_E30N30.TIF"), proxy = TRUE)
tile3040 <- read_stars(file.path(path_dem, "eu_dem_v11_E30N40.TIF"), proxy = TRUE)
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

Visualize DEM data to compare with UK boundaries

The DEM data from Copernicus is in an curvilinear projection, so we need to transform the uk boundaries to that projection in order to adequately filter the data.

uk_sph <- st_transform(uk, crs = st_crs(tile3030))
bbox_sph <- st_bbox(uk_sph)

par(mar = c(2, 2, 2, 0))
plot(st_as_sfc(bbox_sph), graticule = TRUE, axes = TRUE)
plot(tile3030, add = TRUE, main = "Elevation (DEM)")
plot(tile3040, add = TRUE, main = NULL)
plot(uk_sph, add = TRUE, border = 2, lwd = 1.5)
plot(st_as_sfc(bbox_sph), add = TRUE, lwd = 2)

Crop elevation tiles according to UK bounding box

tile3040_sub <- tile3040[bbox_sph]
tile3030_sub <- tile3030[bbox_sph]

tile3040_tmp <- tempfile(fileext = ".tif")
tile3030_tmp <- tempfile(fileext = ".tif")

write_stars(tile3040_sub, tile3040_tmp)

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

write_stars(tile3030_sub, tile3030_tmp)

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

tile3030_sub <- read_stars(tile3030_tmp, proxy = TRUE)
tile3040_sub <- read_stars(tile3040_tmp, proxy = TRUE)

par(mar = c(2, 2, 2, 0))
plot(st_as_sfc(bbox_sph), graticule = TRUE, axes = TRUE)
plot(tile3030_sub, add = TRUE, main = "Elevation (DEM)")
plot(tile3040_sub, add = TRUE, main = NULL)
plot(uk_sph, add = TRUE, border = 2, lwd = 1.5)
plot(st_as_sfc(bbox_sph), add = TRUE, lwd = 2)

Merge tiles of elevation

tiles_merged <- st_mosaic(tile3040_sub, tile3030_sub)
system.time(write_stars(tiles_merged, file.path(path_cleaned, "elev.tif")))

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

#>    user  system elapsed 
#>  87.557  34.034 122.937

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#> 286.690  85.780 375.307