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