tests/testthat/test-postProcess.R

In reproducible: Enhance Reproducibility of R Code

test_that("prepInputs doesn't work (part 3)", {
  skip_on_cran() # too long
  skip_if(getRversion() < "4.1" && isWindows()) # old Windows is failing; not going to fix tests for those
  skip_if_not_installed("sf")
  testInit(c("terra", "sf"),
    tmpFileExt = c(".tif", ".tif", ".tif"),
    opts = list(
      "rasterTmpDir" = tempdir2(rndstr(1, 6)),
      "reproducible.inputPaths" = NULL,
      "reproducible.overwrite" = TRUE,
      "rgdal_show_exportToProj4_warnings" = "none"
    ) # https://gis.stackexchange.com/questions/390945/importing-raster-files-warning-and-extracting-covariates-error-with-raster-and
  )

  options("reproducible.cachePath" = tmpdir)

  # Add a study area to Crop and Mask to
  # Create a "study area"

  coords <- structure(c(-122.98, -116.1, -99.2, -106, -122.98, 59.9, 65.73, 63.58, 54.79, 59.9),
    .Dim = c(5L, 2L)
  )
  coords2 <- structure(c(-115.98, -116.1, -99.2, -106, -122.98, 59.9, 65.73, 63.58, 54.79, 59.9),
    .Dim = c(5L, 2L)
  )

  StudyArea <- terra::vect(coords, "polygons")
  terra::crs(StudyArea) <- crsToUse
  StudyArea2 <- terra::vect(coords2, "polygons")
  terra::crs(StudyArea2) <- crsToUse

  # coords <- structure(c(-122.98, -116.1, -99.2, -106, -122.98, 59.9, 65.73, 63.58, 54.79, 59.9),
  #                     .Dim = c(5L, 2L))
  # coords2 <- structure(c(-115.98, -116.1, -99.2, -106, -122.98, 59.9, 65.73, 63.58, 54.79, 59.9),
  #                      .Dim = c(5L, 2L))
  # Sr1 <- Polygon(coords)
  # Srs1 <- Polygons(list(Sr1), "s1")
  # StudyArea <- SpatialPolygons(list(Srs1), 1L)
  # crs(StudyArea) <- crsToUse
  #
  # Sr1 <- Polygon(coords2)
  # Srs1 <- Polygons(list(Sr1), "s1")
  # StudyArea2 <- SpatialPolygons(list(Srs1), 1L)
  # crs(StudyArea2) <- crsToUse

  nonLatLongProj <- paste(
    "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95",
    "+x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"
  )
  nc <- sf::st_as_sf(StudyArea) # system.file("shape/nc.shp", package="sf"))
  nc1 <- sf::st_transform(nc, nonLatLongProj)
  ncSmall <- sf::st_as_sf(StudyArea2)
  ncSmall <- sf::st_transform(ncSmall, nonLatLongProj)
  ncSmall <- sf::st_buffer(ncSmall, dist = -10000)
  b <- postProcess(nc1, studyArea = ncSmall, filename2 = NULL)
  expect_true(is(b, "sf"))
  expect_equal(terra::ext(b), terra::ext(ncSmall))
  expect_true(sf::st_area(b) < sf::st_area(nc1))

  r <- suppressWarnings(terra::rast(nc1, res = 1000)) # TODO: temporary until raster crs fixes

  rB <- suppressWarnings(terra::rast(nc1, res = 4000)) # TODO: temporary until raster crs fixes
  rSmall <- suppressWarnings(terra::rast(ncSmall, res = 4000)) # TODO: temporary until raster crs fixes

  # Tests with RasterBrick
  r2 <- r1 <- rB
  r1[] <- runif(terra::ncell(rB))
  r2[] <- runif(terra::ncell(rB))

  b <- c(r1, r2)
  terra::crs(b) <- sf::st_crs(ncSmall)$input
  b1 <- postProcess(b, studyArea = ncSmall, useCache = FALSE)
  expect_true(inherits(b1, "SpatRaster"))

  s <- c(r1, r2)
  crs(s) <- crs(nonLatLongProj)
  s1 <- postProcess(s, studyArea = ncSmall, useCache = FALSE)
  expect_true(inherits(s1, "SpatRaster"))
  expect_equal(s1[], b1[], ignore_attr = TRUE)

  b <- writeRaster(b, filename = tmpfile[1], overwrite = TRUE)
  b1 <- postProcess(b, studyArea = ncSmall, useCache = FALSE, filename2 = tmpfile[2], overwrite = TRUE)
  expect_true(inherits(b1, "SpatRaster"))

  s1 <- postProcess(s, studyArea = ncSmall, useCache = FALSE, filename2 = tmpfile[2], overwrite = TRUE)
  expect_true(inherits(s1, "SpatRaster"))

  # Test datatype setting
  dt1 <- "INT2U"
  s <- writeRaster(s, filename = tmpfile[2], overwrite = TRUE)
  s1 <- postProcess(s,
    studyArea = ncSmall, useCache = FALSE, filename2 = tmpfile[1], overwrite = TRUE,
    datatype = dt1
  )
  expect_identical(terra::datatype(s1), rep(dt1, terra::nlyr(s)))

  # Test datatype setting
  dt1 <- c("INT2U", "INT4U")
  s <- writeRaster(s, filename = tmpfile[1], overwrite = TRUE)
  warns <- capture_error({
    s1 <- postProcess(s,
      studyArea = ncSmall, useCache = FALSE, filename2 = tmpfile[2], overwrite = TRUE,
      datatype = dt1
    )
  })
  expect_true(any(grepl("Expecting a single", warns)))

  dt1 <- "INT4U"
  b <- writeRaster(b, filename = tmpfile[2], overwrite = TRUE)
  b1 <- postProcess(b,
    studyArea = ncSmall, useCache = FALSE, filename2 = tmpfile[1], overwrite = TRUE,
    datatype = dt1
  )
  expect_identical(terra::datatype(b1), rep(dt1, terra::nlyr(b1)))

  # now raster with sf
  skip_if_not_installed("terra")
  r1 <- terra::rasterize(terra::vect(nc1), r)
  r2 <- postProcess(r1, studyArea = ncSmall, filename2 = NULL)
  expect_true(is(r2, "SpatRaster"))
  expect_true(terra::ncell(r2) < terra::ncell(r1))
  expect_true((terra::xmin(terra::ext(ncSmall)) - terra::xmin(r2)) < terra::res(r2)[1] * 2)
  expect_true((terra::ymin(terra::ext(ncSmall)) - terra::ymin(r2)) < terra::res(r2)[2] * 2)
  expect_true((terra::ymax(terra::ext(ncSmall)) - terra::ymax(r2)) > -(terra::res(r2)[2] * 2))
  expect_true((terra::xmax(terra::ext(ncSmall)) - terra::xmax(r2)) > -(terra::res(r2)[2] * 2))

  # postProcess
  expect_error(postProcess(1), regexp = "from must be a")
  expect_error(postProcess(list(1, 1)), regexp = "from must be a")

  nc2 <- postProcess(nc1, studyArea = as(ncSmall, "sf"))
  expect_equal(st_area(nc2), st_area(ncSmall))

  # cropInputs
  expect_true(identical(1, cropInputs(1)))
  nonLatLongProj2 <- paste(
    "+proj=lcc +lat_1=51 +lat_2=77 +lat_0=0 +lon_0=-95",
    "+x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"
  )
  nc3 <- suppressWarningsSpecific(
    {
      sf::st_transform(nc1, CRSobj = nonLatLongProj2)
    },
    falseWarnings = "Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition|PROJ support is provided by the sf and terra packages among others"
  )
  nc4 <- cropInputs(nc3, studyArea = ncSmall)
  ncSmall2 <- suppressWarningsSpecific(
    {
      sf::st_transform(ncSmall, CRSobj = nonLatLongProj2)
    },
    falseWarnings = "Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition|PROJ support is provided by the sf and terra packages among others"
  )
  expect_true(isTRUE(all.equal(terra::ext(nc4), terra::ext(ncSmall2))))

  mess <- capture_error({
    nc4 <- cropInputs(nc3, studyArea = 1)
  })
  expect_true(grepl("anyNA.+is not TRUE", mess))

  ncSmallShifted <- ncSmall + 10000000
  ncSmallShifted <- st_as_sf(ncSmallShifted)
  st_crs(ncSmallShifted) <- st_crs(ncSmall)
  mess <- capture_messages(
    cropInputs(ncSmall, studyArea = ncSmallShifted)
  )
  expect_true(any(grepl("overlap", mess)))

  # cropInputs.sf
  nc3 <- st_transform(nc1, crs = nonLatLongProj2)
  nc4 <- cropInputs(nc3, studyArea = ncSmall)
  ncSmall2 <- st_transform(ncSmall, crs = nonLatLongProj2)
  expect_true(isTRUE(all.equal(terra::ext(nc4), terra::ext(ncSmall2))))

  # studyArea as spatial object
  nc5 <- cropInputs(nc3, studyArea = ncSmall)
  ncSmall2 <- st_transform(ncSmall, crs = nonLatLongProj2)
  expect_true(isTRUE(all.equal(terra::ext(nc5), terra::ext(ncSmall2))))
  expect_true(isTRUE(all.equal(terra::ext(nc5), terra::ext(nc4))))


  # rasterToMatch
  nc5 <- cropTo(nc3, cropTo = r)
  nc5Extent_r <- st_transform(nc5, crs = crs(r))
  expect_true(isTRUE(abs(terra::xmin(r) - sf::st_bbox(nc5Extent_r)["xmin"]) < terra::res(r)[1]))
  expect_true(isTRUE(abs(terra::ymin(r) - sf::st_bbox(nc5Extent_r)["ymin"]) < terra::res(r)[1]))
  expect_true(isTRUE(abs(terra::xmax(r) - sf::st_bbox(nc5Extent_r)["xmax"]) < terra::res(r)[1]))
  expect_true(isTRUE(abs(terra::ymax(r) - sf::st_bbox(nc5Extent_r)["ymax"]) < terra::res(r)[1]))

  ncSmallShifted <- ncSmall + 10000000
  ncSmallShifted <- st_as_sf(ncSmallShifted)
  st_crs(ncSmallShifted) <- st_crs(ncSmall)
  expect_message(
    regexp = "results in no data",
    aaa <- cropInputs(ncSmall, studyArea = ncSmallShifted)
  )
  expect_s3_class(aaa, "sf")
  # expect_true(NROW(out11) == 0)

  # LINEARRING Example
  p6 <- terra::vect("POLYGON ((0 60, 0 0, 60 0, 60 20, 100 20, 60 20, 60 60, 0 60))")
  p6a <- fixErrorsIn(p6)
  expect_true(terra::is.valid(p6a))
  expect_false(terra::is.valid(p6))
  # projectInputs pass through
  expect_error(projectInputs(x = 1), "argument .+ is missing")
})

test_that("writeOutputs with non-matching filename2", {
  testInit(c("terra"), tmpFileExt = c(".grd", ".tif"))

  r <- terra::rast(terra::ext(0, 10, 0, 10), vals = rnorm(100))
  r <- terra::writeRaster(r, file = tmpfile[1], overwrite = TRUE)
  r[] <- r[]
  warn <- capture_warnings({
    r1 <- writeOutputs(r, filename2 = tmpfile[2])
  })
  r2 <- terra::rast(Filenames(r1))
  vals1 <- r2[]
  vals2 <- r1[]
  vals3 <- r[]
  expect_true(identical(vals1, vals2))
  expect_true(identical(vals1, vals3))
  expect_false(identical(normPath(Filenames(r)), normPath(Filenames(r1))))
  expect_true(identical(normPath(Filenames(r2)), normPath(Filenames(r1))))
})

test_that("cropInputs crops too closely when input projections are different", {
  skip_on_cran()
  # skip_on_ci() ## TODO: why is this failing on GHA but not locally??? (2022-11-04)

  testInit("terra", opts = list(
    "rasterTmpDir" = tempdir2(rndstr(1, 6)),
    "reproducible.overwrite" = TRUE,
    "reproducible.inputPaths" = NULL,
    "rgdal_show_exportToProj4_warnings" = "none" # https://gis.stackexchange.com/questions/390945/importing-raster-files-warning-and-extracting-covariates-error-with-raster-and
  ), needGoogleDriveAuth = TRUE)

  ext2 <- terra::ext(c(
    xmin = -3229772.32501426,
    xmax = 3680227.67498574,
    ymin = -365833.605586135,
    ymax = 3454166.39441387
  ))
  x <- terra::rast(ext2,
    crs = paste(
      "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0",
      "+a=6370997 +b=6370997 +units=m +no_defs"
    ),
    res = c(10000, 10000)
  )
  x <- terra::setValues(x, 1)

  RTMext <- terra::ext(c(
    xmin = -1613500.00000023,
    xmax = -882500.000000228,
    ymin = 7904500,
    ymax = 9236000
  ))
  RTM <- terra::rast(RTMext,
    crs = paste(
      "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95 +x_0=0 +y_0=0",
      "+ellps=GRS80 +units=m +no_defs"
    ),
    res = c(250, 250)
  )
  RTM <- setValues(RTM, 2)
  out <- postProcess(x = x, rasterToMatch = RTM, filename2 = NULL)
  expect_true(nrow(out[is.na(out[]) & is.na(RTM[])]) == 0)
})

test_that("maskInputs errors when x is Lat-Long", {
  skip_on_cran()
  skip_on_ci()
  testInit("sf", opts = list(
    "rasterTmpDir" = tempdir2(rndstr(1, 6)),
    "reproducible.overwrite" = TRUE,
    "reproducible.inputPaths" = NULL
  ), needGoogleDriveAuth = TRUE)

  i <- 0
  roads <- list()
  clearCache()

  smallSA <- terra::ext(c(
    -117.580383168455, -117.43120279669,
    61.0576330401172, 61.0937107807574
  ))

  x <- terra::rast(smallSA,
    crs = "+proj=longlat +ellps=GRS80 +no_defs",
    res = c(0.001, 0.001)
  )
  suppressWarnings(smallSA <- terra::vect(terra::ext(x), "polygons"))
  terra::crs(smallSA) <- terra::crs(x)

  noisyOutput <- capture.output(
    suppressWarningsSpecific(
      falseWarnings = "attribute variables|st_buffer does not correctly buffer longitude",
      roads1 <- prepInputs(
        targetFile = "miniRoad.shp",
        alsoExtract = "similar",
        url = "https://drive.google.com/file/d/1Z6ueq8yXtUPuPWoUcC7_l2p0_Uem34CC",
        studyArea = smallSA,
        useCache = FALSE,
        fun = "sf::st_read",
        destinationPath = tmpdir,
        filename2 = "miniRoad.shp"
      )
    )
  )
  # clearCache()
  noisyOutput <- capture.output(
    roads2 <- prepInputs(
      targetFile = "miniRoad.shp",
      alsoExtract = "similar",
      url = "https://drive.google.com/file/d/1Z6ueq8yXtUPuPWoUcC7_l2p0_Uem34CC",
      # studyArea = smallSA,
      useCache = FALSE,
      fun = "sf::st_read",
      destinationPath = tmpdir,
      filename2 = "miniRoads"
    )
  )
  attr(roads1, "tags") <- NULL

  # There are floating point issues with 32 bit vs 64 bit approaches. The following fails:
  # expect_true(all.equal(roads[[1]], roads[[2]], check.attributes = FALSE))

  # diffs <- sum(abs(unlist(lapply(sf::st_geometry(roads[[1]]), as.numeric)) -
  #                    unlist(lapply(sf::st_geometry(roads[[2]]), as.numeric))))
  # expect_true(diffs < 0.0001)
  expect_true(terra::compareGeom(terra::rast(terra::ext(roads1)), terra::rast(terra::ext(smallSA))))
  expect_error(terra::compareGeom(terra::rast(terra::ext(roads2)), terra::rast(terra::ext(smallSA))))
  expect_true(terra::ext(roads2) > terra::ext(roads1))
})

test_that("prepInputs doesn't work (part 3)", {
  skip("The Google Drive url is dead")
  if (interactive()) {
    testInit(needGoogleDriveAuth = TRUE)

    # Tati's reprex
    # tmpdir <- "/mnt/d/temp/Cache"
    wd <- checkPath(file.path(tmpdir, "reprex"), create = TRUE)
    ranges <- prepInputs(
      url = "https://drive.google.com/file/d/1AfGfRjaDsdq3JqcsidGRo3N66OUjRJnn",
      destinationPath = wd,
      fun = "terra::vect"
    )
    LCC05 <- prepInputs(
      url = "https://drive.google.com/file/d/1g9jr0VrQxqxGjZ4ckF6ZkSMP-zuYzHQC",
      targetFile = "LCC2005_V1_4a.tif",
      studyArea = ranges,
      fun = "terra::rast",
      destinationPath = wd
    )
    sumNonNAs <- sum(!is.na(!LCC05[]))

    # These are suitably vague that they will capture the mask if it gets it right
    expect_true(sumNonNAs < 38000000)
    expect_true(sumNonNAs > 37000000)
  }
})