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, resolution = 1000)) # TODO: temporary until raster crs fixes
rB <- suppressWarnings(terra::rast(nc1, resolution = 4000)) # TODO: temporary until raster crs fixes
rSmall <- suppressWarnings(terra::rast(ncSmall, resolution = 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, to = r2), regexp = "from must be a")
expect_error(postProcess(list(1, 1), to = r2), 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, filename = 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"
),
resolution = 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"
),
resolution = 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",
resolution = 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)
}
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.