# \donttest{
## Set up file paths
td <- tempdir()
dem_file <- file.path(td, "dem.csv")
vrt_header_file <- file.path(td, "tmp.vrt")
out_raster <- file.path(td, "tmp.tiff")
## Create file of points with x-, y-, and z-coordinates
pts <-
data.frame(Easting = c(86943.4, 87124.3, 86962.4, 87077.6),
Northing = c(891957, 892075, 892321, 891995),
Elevation = c(139.13, 135.01, 182.04, 135.01))
write.csv(pts, file = dem_file, row.names = FALSE)
## Prepare a matching VRT file
vrt_header <- c(
'',
' ',
paste0(' ',dem_file,''),
' wkbPoint',
' ',
' ',
''
)
cat(vrt_header, file = vrt_header_file, sep = "\n")
## Test it out
gdal_grid(src_datasource = vrt_header_file,
dst_filename = out_raster,
a = "invdist:power=2.0:smoothing=1.0",
txe = c(85000, 89000), tye = c(894000, 890000),
outsize = c(400, 400),
of = "GTiff", ot = "Float64", l = "dem")
## Check that it works
if(requireNamespace("terra", quietly = TRUE)) {
library(terra)
plot(rast(out_raster))
text(Northing ~ Easting, data = pts,
labels = seq_len(nrow(pts)), cex = 0.7)
}
# }
Run the code above in your browser using DataLab