Nothing
test_that(".getGDALformat works", {
expect_equal(.getGDALformat("test.tif"), "GTiff")
expect_equal(.getGDALformat("test.TIF"), "GTiff")
expect_equal(.getGDALformat("test.img"), "HFA")
expect_equal(.getGDALformat("test.vrt"), "VRT")
expect_equal(.getGDALformat("test.VRT"), "VRT")
expect_null(.getGDALformat("test.unknown"))
})
test_that("calc writes correct results", {
# bioclimatic index
elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
expr <- "round( ((ELEV_M * 3.281 - 5449) / 100) +
((pixelLat - 42.16) * 4) +
((-116.39 - pixelLon) * 1.25) )"
hi_file <- calc(expr = expr,
rasterfiles = elev_file,
var.names = "ELEV_M",
dtName = "Int16",
setRasterNodataValue = TRUE)
on.exit(deleteDataset(hi_file))
ds <- new(GDALRaster, hi_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 2978)
# NDVI
b4_file <- system.file("extdata/sr_b4_20200829.tif", package="gdalraster")
b5_file <- system.file("extdata/sr_b5_20200829.tif", package="gdalraster")
expr <- "((B5-B4)/(B5+B4)) * 1000" # rescale for checksum
ndvi_file <- calc(expr = expr,
rasterfiles = c(b4_file, b5_file),
var.names = c("B4", "B5"),
dtName = "Float32",
nodata_value = -32767,
setRasterNodataValue = TRUE)
on.exit(deleteDataset(ndvi_file))
ds <- new(GDALRaster, ndvi_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 63632)
# recode
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
tif_file <- paste0(tempdir(), "/", "storml_lndscp.tif")
on.exit(deleteDataset(tif_file))
createCopy("GTiff", tif_file, lcp_file)
expr <- "ifelse( SLP >= 40 & FBFM %in% c(101,102), 99, FBFM)"
calc(expr = expr,
rasterfiles = c(lcp_file, lcp_file),
bands = c(2, 4),
var.names = c("SLP", "FBFM"),
dstfile = tif_file,
out_band = 4,
write_mode = "update")
ds <- new(GDALRaster, tif_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 28017)
# multiband output
# https://github.com/USDAForestService/gdalraster/issues/319
expr <- "
to_terrainrgb <- function(dtm) {
startingvalue <- 10000
precision <- 0.1
rfactor <- 256*256 * precision
gfactor <- 256 * precision
r <- floor((startingvalue +dtm)*(1/precision) / 256 / 256)
g <- floor((startingvalue +dtm - r*rfactor)*(1/precision) / 256)
b <- floor((startingvalue +dtm - r*rfactor - g*gfactor)*(1/precision))
return(c(r,g,b))
}
to_terrainrgb(ELEV)"
out_file <- "/vsimem/multiband-calc.tif"
result <- calc(expr = expr,
rasterfiles = elev_file,
var.names = "ELEV",
dstfile = out_file,
dtName = "Byte",
out_band = 1:3,
nodata_value = 0)
expect_equal(result, out_file)
# revert out_file using from_terrainrgb()
expr <- "
from_terrainrgb <- function(R,G,B) {
elevation <- -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
return(elevation)
}
from_terrainrgb(R,G,B)"
revert_file <- "/vsimem/multiband-calc-revert.tif"
result <- calc(expr = expr,
rasterfiles = c(out_file, out_file, out_file),
bands = c(1, 2 ,3),
var.names = c("R", "G", "B"),
dstfile = revert_file,
dtName = "Int16",
nodata_value = -10000,
setRasterNodataValue = TRUE)
expect_equal(result, revert_file)
# compare revert_file with the original elev_file
ds <- new(GDALRaster, elev_file)
orig_elev_values <- read_ds(ds)
attributes(orig_elev_values) <- NULL
ds$close()
ds2 <- new(GDALRaster, revert_file)
revert_elev_values <- read_ds(ds2)
attributes(revert_elev_values) <- NULL
ds2$close()
expect_equal(revert_elev_values, orig_elev_values)
expect_equal(mean(revert_elev_values, na.rm = TRUE),
mean(orig_elev_values, na.rm = TRUE))
deleteDataset(out_file)
deleteDataset(revert_file)
# test errors from input validation
expr <- "((B5-B4)/(B5+B4))"
expect_error(calc(expr = expr,
rasterfiles = c(b4_file, paste0(b5_file, ".error")),
var.names = c("B4", "B5"),
dtName = "Float32",
nodata_value = -32767,
setRasterNodataValue = TRUE))
expect_error(calc(expr = expr,
rasterfiles = c(b4_file, b5_file),
bands = 1,
var.names = c("B4", "B5"),
dtName = "Float32",
nodata_value = -32767,
setRasterNodataValue = TRUE))
expect_error(calc(expr = expr,
rasterfiles = c(b4_file, b5_file),
var.names = c("B4", "B5", "err"),
dtName = "Float32",
nodata_value = -32767,
setRasterNodataValue = TRUE))
expect_error(calc(expr = expr,
rasterfiles = c(b4_file, b5_file),
dtName = "Float32",
nodata_value = -32767,
setRasterNodataValue = TRUE))
expect_error(calc(expr = expr,
rasterfiles = c(b4_file, b5_file),
var.names = c("B4", "B5"),
dstfile = "unknown_file_type.err",
dtName = "Float32",
nodata_value = -32767,
setRasterNodataValue = TRUE))
expect_error(calc(expr = expr,
rasterfiles = c(b4_file, b5_file),
var.names = c("B4", "B5"),
dtName = "Int64",
setRasterNodataValue = TRUE))
expect_error(calc(expr = expr,
rasterfiles = c(b4_file, elev_file),
var.names = c("B4", "B5"),
dtName = "Float32",
nodata_value = -32767,
setRasterNodataValue = TRUE))
})
test_that("combine writes correct output", {
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
rasterfiles <- c(lcp_file, lcp_file)
bands <- c(4, 5)
var.names <- c("fbfm", "tree_cov")
cmb_file <- paste0(tempdir(), "/", "fbfm_cov_cmbid.tif")
on.exit(deleteDataset(cmb_file))
df <- combine(rasterfiles, var.names, bands, cmb_file)
expect_equal(nrow(df), 29)
df <- NULL
ds <- new(GDALRaster, cmb_file, TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 43024)
evt_file <- system.file("extdata/storml_evt.tif", package="gdalraster")
df <- combine(evt_file)
expect_equal(nrow(df), 24)
})
test_that("rasterFromRaster works", {
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
slpp_file <- paste0(tempdir(), "/", "storml_slpp.tif")
on.exit(deleteDataset(slpp_file))
options = c("COMPRESS=LZW")
rasterFromRaster(srcfile = lcp_file,
dstfile = slpp_file,
nbands = 1,
dtName = "Int16",
options = options,
init = -32767)
ds <- new(GDALRaster, slpp_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 47771)
slpp_file_img <- paste0(tempdir(), "/", "storml_slpp.img")
on.exit(deleteDataset(slpp_file_img))
options = c("COMPRESSED=YES")
rasterFromRaster(srcfile = lcp_file,
dstfile = slpp_file_img,
nbands = 1,
dtName = "Int16",
options = options,
init = -32767)
ds <- new(GDALRaster, slpp_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 47771)
slpp_file_envi <- paste0(tempdir(), "/", "storml_slpp")
on.exit(deleteDataset(slpp_file_envi))
# error on unknown file format:
expect_error(rasterFromRaster(srcfile = lcp_file,
dstfile = slpp_file_envi,
nbands = 1,
dtName = "Int16",
init = -32767))
rasterFromRaster(srcfile = lcp_file,
dstfile = slpp_file_envi,
fmt = "ENVI",
nbands = 1,
dtName = "Int16",
init = -32767)
ds <- new(GDALRaster, slpp_file_envi, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 47771)
})
test_that("rasterToVRT works", {
## resample
evt_file <- system.file("extdata/storml_evt.tif", package="gdalraster")
vrt_file <- rasterToVRT(evt_file,
resolution=c(90,90),
resampling="mode")
on.exit(deleteDataset(vrt_file))
ds <- new(GDALRaster, vrt_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 23684)
## clip raster to polygon extent
bnd = "POLYGON ((324467.3 5104814.2, 323909.4 5104365.4, 323794.2
5103455.8, 324970.7 5102885.8, 326420.0 5103595.3, 326389.6 5104747.5,
325298.1 5104929.4, 325298.1 5104929.4, 324467.3 5104814.2))"
# src_align = TRUE
vrt_file <- rasterToVRT(evt_file,
subwindow = bbox_from_wkt(bnd),
src_align=TRUE)
on.exit(deleteDataset(vrt_file))
ds <- new(GDALRaster, vrt_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 18148)
# src_align = FALSE
vrt_file <- rasterToVRT(evt_file,
subwindow = bbox_from_wkt(bnd),
src_align=FALSE)
on.exit(deleteDataset(vrt_file))
ds <- new(GDALRaster, vrt_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 17631)
# subwindow outside raster extent
expect_error(rasterToVRT(evt_file,
subwindow = bbox_from_wkt(g_buffer(bnd, 10000)),
src_align=TRUE))
## subset and pixel align two rasters
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
ds_lcp <- new(GDALRaster, lcp_file, read_only=TRUE)
b5_file <- system.file("extdata/sr_b5_20200829.tif", package="gdalraster")
vrt_file <- rasterToVRT(b5_file,
resolution = ds_lcp$res(),
subwindow = ds_lcp$bbox(),
src_align = FALSE)
on.exit(deleteDataset(vrt_file))
ds_lcp$close()
ds <- new(GDALRaster, vrt_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 49589)
## kernel filter
elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
krnl <- c(0.11111, 0.11111, 0.11111,
0.11111, 0.11111, 0.11111,
0.11111, 0.11111, 0.11111)
vrt_file <- rasterToVRT(elev_file, krnl=krnl)
on.exit(deleteDataset(vrt_file))
ds <- new(GDALRaster, vrt_file, read_only=TRUE)
dm <- ds$dim()
chk <- ds$getChecksum(1, 0, 0, dm[1], dm[2])
ds$close()
expect_equal(chk, 46590)
expect_error(rasterToVRT(elev_file, resolution=c(90,90), krnl=krnl))
})
test_that("dem_proc runs without error", {
elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
hs_file <- paste0(tempdir(), "/", "storml_hillshade.tif")
on.exit(deleteDataset(hs_file))
expect_true(dem_proc("hillshade", elev_file, hs_file))
sq <- seq(2400, 3100, by = 100)
col <- col2rgb(hcl.colors(length(sq)))
pal <- tempfile(fileext = ".txt")
on.exit(unlink(pal))
writeLines(paste(sq, col[1, ], col[2, ], col[3, ]), pal)
cr_file <- tempfile(fileext = ".tif")
on.exit(deleteDataset(cr_file))
expect_true(dem_proc("color-relief", elev_file, cr_file, color_file = pal))
})
test_that("polygonize runs without error", {
evt_file <- system.file("extdata/storml_evt.tif", package="gdalraster")
dsn <- paste0(tempdir(), "/", "storml_evt.shp")
layer <- "storml_evt"
fld <- "evt_value"
expect_true(polygonize(evt_file, dsn, layer, fld))
# overwrite
expect_true(polygonize(evt_file, dsn, layer, fld, connectedness=8,
overwrite=TRUE))
# with mask_file
expr <- "ifelse(EVT == -9999, 0, 1)"
mask_file <- calc(expr, rasterfiles=evt_file, var.names="EVT")
expect_true(polygonize(evt_file, dsn, layer, fld, mask_file=mask_file,
overwrite=TRUE))
# nomask
expect_true(polygonize(evt_file, dsn, layer, fld, nomask=TRUE,
overwrite=TRUE))
# append: layer already exists so expect warning if lco given
expect_warning(polygonize(evt_file, dsn, layer, fld, connectedness=8,
lco="SPATIAL_INDEX=YES"))
deleteDataset(dsn)
deleteDataset(mask_file)
# GPKG
set_config_option("SQLITE_USE_OGR_VFS", "YES")
set_config_option("OGR_SQLITE_JOURNAL", "MEMORY")
dsn <- paste0(tempdir(), "/", "storml_evt.gpkg")
layer <- "lf_evt"
opt <- "VERSION=1.3"
expect_true(polygonize(evt_file, dsn, layer, fld, dsco=opt))
# existing out_dsn, but create a new layer
layer <- "lf_evt_8connect"
opt <- c("GEOMETRY_NULLABLE=NO","DESCRIPTION=LF EVT 8-connected polygons")
expect_true(polygonize(evt_file,dsn,layer,fld,connectedness=8,lco=opt))
set_config_option("SQLITE_USE_OGR_VFS", "")
set_config_option("OGR_SQLITE_JOURNAL", "")
deleteDataset(dsn)
# unrecognized output format
dsn <- paste0(tempdir(), "/", "storml_evt.txt")
expect_error(polygonize(evt_file, dsn, layer, fld))
})
test_that("rasterize runs without error", {
# layer from sql query
dsn <- system.file("extdata/ynp_fires_1984_2022.gpkg", package="gdalraster")
sql <- "SELECT * FROM mtbs_perims ORDER BY mtbs_perims.ig_year"
out_file <- paste0(tempdir(), "/", "ynp_fires_1984_2022.tif")
res <- rasterize(src_dsn = dsn,
dstfile = out_file,
sql = sql,
burn_attr = "ig_year",
tr = c(90,90),
tap = TRUE,
dtName = "Int16",
dstnodata = -9999,
init = -9999,
co = c("TILED=YES","COMPRESS=LZW"))
expect_true(res)
ds <- new(GDALRaster, out_file, TRUE)
bbox <- ds$bbox()
dm <- ds$dim()
ds$close()
deleteDataset(out_file)
# layer with where clause
out_file <- paste0(tempdir(), "/", "ynp_fires_1988.tif")
layer <- "mtbs_perims"
where <- "ig_year = 1988"
res <- rasterize(src_dsn = dsn,
dstfile = out_file,
layer = layer,
where = where,
burn_value = 1,
te = bbox,
ts = dm[1:2],
dtName = "Byte",
init = 0,
fmt = "GTiff",
co = c("TILED=YES","COMPRESS=LZW"),
add_options = "-q")
expect_true(res)
deleteDataset(out_file)
# update existing raster in-place
# create the vector layer
dsn <- tempfile(fileext = ".geojson")
lyr <- ogr_ds_create("GeoJSON", dsn, "OGRGeoJSON", geom_type = "POLYGON",
return_obj = TRUE)
pts <- matrix(c(0.25, 0.25, 0.75, 0.25, 0.75, 0.75, 0.25, 0.75, 0.25, 0.25),
ncol = 2, byrow = TRUE)
feat <- list()
feat$geometry <- g_create("POLYGON", pts)
lyr$createFeature(feat)
lyr$close()
# create the destination raster
ds <- create("GTiff", "/vsimem/test.tif", xsize = 100, ysize = 100,
nbands = 1, dataType = "Byte", options = "COMPRESS=DEFLATE",
return_obj = TRUE)
ds$setGeoTransform(c(0.0, 0.01, 0.0, 1.0, 0.0, -0.01))
ds$setProjection(epsg_to_wkt(4326))
ds$fillRaster(1, 0, 0)
# update the destination pixels in-place
res <- rasterize(dsn, ds, burn_value = 1)
expect_true(res)
r <- read_ds(ds)
expect_equal(sum(r), 2500)
unlink(dsn)
ds$close()
vsi_unlink("/vsimem/test.tif")
})
test_that("pixel_extract wrapper returns correct data", {
# the C++ class method GDALRaster::pixel_extract() is tested in
# test-GDALRaster-class.R
pt_file <- system.file("extdata/storml_pts.csv", package="gdalraster")
pts <- read.csv(pt_file)
raster_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
# single pixel extract
extr <- pixel_extract(raster_file, pts[-1])
colnames(extr) <- NULL
dim(extr) <- NULL
expected_values <- c(2648, 2865, 2717, 2560, 2916,
2633, 2548, 2801, 2475, 2822)
expect_equal(extr, expected_values)
# as a GDALRaster object
ds <- new(GDALRaster, raster_file)
extr_ds <- pixel_extract(ds, pts[-1])
colnames(extr_ds) <- NULL
dim(extr_ds) <- NULL
expect_equal(extr_ds, expected_values)
# with missing values in the input
pts_na <- rbind(pts, c(11, NA_real_, NA_real_))
extr_na <- pixel_extract(ds, pts_na[-1])
colnames(extr_na) <- NULL
dim(extr_na) <- NULL
expect_na <- c(expected_values, NA_real_)
expect_equal(extr_na, expect_na)
# with NaN in the input
pts_nan <- rbind(pts, c(11, NaN, NaN))
extr_nan <- pixel_extract(ds, pts_na[-1])
colnames(extr_nan) <- NULL
dim(extr_nan) <- NULL
expect_nan <- c(expected_values, NA_real_)
expect_equal(extr_nan, expect_nan)
ds$close()
# interpolated values
extr_bilinear <- pixel_extract(raster_file, pts[-1], interp = "bilinear")
colnames(extr_bilinear) <- NULL
dim(extr_bilinear) <- NULL
expected_values <- c(2649.217, 2881.799, 2716.290, 2558.797, 2920.404,
2629.495, 2548.250, 2810.543, 2478.609, 2819.776)
expect_equal(extr_bilinear, expected_values, tolerance = 1e-4)
# kernel values
extr_3x3 <- pixel_extract(raster_file, pts[-1], krnl_dim = 3)
colnames(extr_3x3) <- NULL
dimnames(extr_3x3) <- NULL
expected_values <- c(
2660, 2654, 2651, 2652, 2648, 2646, 2653, 2647, 2645,
2924, 2895, 2901, 2893, 2865, 2869, 2863, 2838, 2841,
NA_real_, 2709, 2704, NA_real_, 2717, 2717, NA_real_, 2725, 2728,
2554, 2558, 2562, 2557, 2560, 2564, 2559, 2562, 2566,
2932, 2920, 2908, 2927, 2916, 2903, 2923, 2911, 2899,
2633, 2642, 2650, 2627, 2633, 2640, 2619, 2624, 2630,
2548, 2548, 2550, 2549, 2548, 2550, 2550, 2548, 2550,
2785, 2777, 2776, 2807, 2801, 2803, 2833, 2825, 2827,
2489, 2482, 2475, 2479, 2475, 2470, 2474, 2471, 2468,
2832, 2810, 2789, 2839, 2822, 2800, 2833, 2811, 2788)
expected_values <- matrix(expected_values, nrow = 10, ncol = 9,
byrow = TRUE)
expect_equal(extr_3x3, expected_values)
# transform the xy
ds$open(read_only = TRUE)
pts_nad83 <- transform_xy(pts[-1], ds$getProjection(), "NAD83")
extr <- pixel_extract(raster_file, pts_nad83, xy_srs = "NAD83")
colnames(extr) <- NULL
dim(extr) <- NULL
expected_values <- c(2648, 2865, 2717, 2560, 2916,
2633, 2548, 2801, 2475, 2822)
expect_equal(extr, expected_values)
ds$close()
# input validation
expect_error(pixel_extract(ds, pts[-1])) # ds closed
expect_error(pixel_extract(NULL, pts[-1]))
expect_error(pixel_extract(raster_file))
expect_error(pixel_extract(raster_file, as.character(pts[-1])))
expect_error(pixel_extract(raster_file, pts)) # more than two columns
pts[, 2] <- as.character(pts[, 1])
expect_error(pixel_extract(raster_file, pts[-1]))
expect_error(pixel_extract(raster_file, pts[-1], bands = "1"))
expect_error(pixel_extract(raster_file, pts[-1], interp = 2))
expect_error(pixel_extract(raster_file, pts[-1], krnl_dim = c(1, 1, 1)))
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.