Data processing: testing warping methods on elevation

We have different datasets coming with different projections. Converting one raster to another grid and projection is called warping. Here we test different warping methods to check the implications (e.g. missing data).

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

source(file.path(path_proj, "src", "stars-aux.R"))

elev <- read_stars(file.path(path_cleaned, "elev.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
#> projected CRS:  OSGB 1936 / British National Grid

elev <- st_dim_reset(st_as_stars(elev[, 5000:7000, 15500:16500]))

uk_bbox_1km <- st_as_stars(st_bbox(elev), dx = 1000, dy = 1000, values = 1)

plot(uk_bbox_1km, reset = FALSE, main = NA)
plot(elev, add = TRUE, border = 1:4)

Visualize current elevation dataset

plot(elev, reset = FALSE, main = NA)
plot(st_transform(uk, st_crs(elev)), add = TRUE, border = 2, lwd = 2)

Change resolution in a new CRS

Method: near

system.time(elev_1km <- st_warp(elev, uk_bbox_1km, use_gdal = TRUE, method =
                                "near"))

#>    user  system elapsed 
#>   0.307   0.032   0.340

plot(elev, reset = FALSE, main = NA, col = hcl.colors(10))
plot(elev_1km, add = TRUE, reset = FALSE, main = NA)
plot(st_transform(uk, st_crs(elev_1km)), add = TRUE, border = 2, lwd = 2)

Method: average

system.time(elev_1km <- st_warp(elev, uk_bbox_1km, use_gdal = TRUE, method =
                                "average"))

#>    user  system elapsed 
#>   0.357   0.036   0.393

plot(elev, reset = FALSE, main = NA, col = hcl.colors(10))
plot(elev_1km, add = TRUE, reset = FALSE, main = NA)
plot(st_transform(uk, st_crs(elev_1km)), add = TRUE, border = 2, lwd = 2)

Method: median

system.time(elev_1km <- st_warp(elev, uk_bbox_1km, use_gdal = TRUE, method =
                                "med"))

#>    user  system elapsed 
#>   0.454   0.084   0.539

plot(elev, reset = FALSE, main = NA, col = hcl.colors(10))
plot(elev_1km, add = TRUE, reset = FALSE, main = NA)
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 
#>  26.020   1.335  27.320