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