tests/testthat/test-gdalraster_proc.R

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

Try the gdalraster package in your browser

Any scripts or data that you put into this service are public.

gdalraster documentation built on June 8, 2025, 12:37 p.m.