test_that("testing terra", { # if (interactive()) { testInit("terra", needGoogleDriveAuth = FALSE, opts = list( reproducible.useMemoise = FALSE, "rgdal_show_exportToProj4_warnings" = "none" ) ) opts <- options(reproducible.cachePath = tmpCache) on.exit( { options(opts) }, add = TRUE ) skip_if_not_installed("terra") f <- system.file("ex/elev.tif", package = "terra") tf <- tempfile(fileext = ".tif") tf1 <- tempfile(fileext = ".tif") tf2 <- tempfile(fileext = ".tif") tf3 <- tempfile(fileext = ".tif") tf4 <- tempfile(fileext = ".tif") tf5 <- tempfile(fileext = ".tif") tf6 <- tempfile(fileext = ".tif") file.copy(f, tf) file.copy(f, tf1) file.copy(f, tf2) file.copy(f, tf3) file.copy(f, tf4) file.copy(f, tf5) file.copy(f, tf6) r <- list(terra::rast(f), terra::rast(tf)) r1 <- list(terra::rast(tf1), terra::rast(tf2)) r2 <- list(terra::rast(tf3), terra::rast(tf4)) rmem <- r2 terra::values(rmem[[1]]) <- values2(rmem[[1]]) terra::values(rmem[[2]]) <- values2(rmem[[2]]) fn <- function(listOf) { listOf } # Test Cache of various nested and non nested SpatRaster # double nest b <- Cache(fn, list(r, r1), cacheRepo = tmpCache) expect_true(is(b, "list")) expect_true(is(b[[1]], "list")) expect_true(is(b[[1]][[1]], "SpatRaster")) # Single nest b <- Cache(fn, r, cacheRepo = tmpCache) expect_true(is(b, "list")) expect_true(is(b[[1]], "SpatRaster")) # mixed nest b <- Cache(fn, list(r[[1]], r1), cacheRepo = tmpCache) expect_true(is(b, "list")) expect_true(is(b[[1]], "SpatRaster")) expect_true(is(b[[2]][[1]], "SpatRaster")) # mix memory and disk b <- Cache(fn, list(r[[1]], r1, rmem), cacheRepo = tmpCache) expect_true(is(b, "list")) expect_true(is(b[[1]], "SpatRaster")) expect_true(is(b[[2]][[1]], "SpatRaster")) expect_true(terra::inMemory(b[[3]][[1]])) expect_true(!terra::inMemory(b[[2]][[1]])) expect_true(!terra::inMemory(b[[1]])) f <- system.file("ex/lux.shp", package = "terra") vOrig <- terra::vect(f) v <- vOrig[1:2, ] rf <- system.file("ex/elev.tif", package = "terra") xOrig <- terra::rast(rf) elevRas <- terra::deepcopy(xOrig) xCut <- terra::classify(xOrig, rcl = 5) xVect <- terra::as.polygons(xCut) xVect2 <- terra::deepcopy(xVect) y <- terra::deepcopy(elevRas) y[y > 200 & y < 300] <- NA_integer_ terra::values(elevRas) <- rep(1L, terra::ncell(y)) vRast <- terra::rast(v, resolution = 0.008333333) # SR, SR t1 <- postProcessTo(elevRas, y) expect_true(sum(is.na(t1[]) != is.na(y[])) == 0) t7 <- postProcessTo(elevRas, projectTo = y) expect_true(identical(t7, elevRas)) t8 <- postProcessTo(elevRas, maskTo = y) expect_true(all.equal(t8, t1)) t9 <- postProcessTo(elevRas, cropTo = vRast) expect_true(terra::ext(v) <= terra::ext(t9)) # SR, SV t2 <- postProcessTo(elevRas, v) # No crop t3 <- postProcessTo(elevRas, maskTo = v) expect_true(terra::ext(t3) == terra::ext(elevRas)) if (.requireNamespace("sf")) { vsf <- sf::st_as_sf(v) vOrigsf <- sf::st_as_sf(vOrig) } t4 <- postProcessTo(elevRas, cropTo = v, maskTo = v) expect_true(terra::ext(t4) == terra::ext(t2)) t5 <- postProcessTo(elevRas, cropTo = v, maskTo = v, projectTo = v) expect_true(identical(t5[], t2[])) if (.requireNamespace("sf")) { t5sf <- postProcessTo(elevRas, cropTo = vsf, maskTo = vsf, projectTo = vsf) expect_true(identical(t5sf[], t2[])) } t6 <- terra::extract(elevRas, v, mean, na.rm = TRUE) expect_true(all(t6$elevation == 1)) expect_true(NROW(t6) == 2) # t10 <- postProcessTo(xVect, v) expect_true(terra::ext(t10) < terra::ext(xVect)) # ## following #253 # https://github.com/PredictiveEcology/reproducible/issues/253#issuecomment-1263562631 tf1 <- tempfile(fileext = ".shp") t11 <- suppressWarnings({ postProcessTo(xVect, v, writeTo = tf1) }) ## WARNING: Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition tw_t11 <- terra::wrap(t11) vv <- terra::vect(tf1) tw_vv <- terra::wrap(vv) expect_true(terra::same.crs(tw_vv@crs, tw_t11@crs)) ## following #253 with different driver ## https://github.com/PredictiveEcology/reproducible/issues/253#issuecomment-1263562631 tf1 <- tempfile(fileext = ".gpkg") t11 <- suppressWarnings({ postProcessTo(xVect, v, writeTo = tf1) }) ## WARNING: GDAL Message 6: dataset does not support layer creation option ENCODING tw_t11 <- terra::wrap(t11) vv <- terra::vect(tf1) tw_vv <- terra::wrap(vv) expect_equivalent(tw_vv, tw_t11) ## TODO: not identical # Test fixErrorTerra v1 <- terra::simplifyGeom(v) gv1 <- terra::geom(v1) gv1[gv1[, "geom"] == 2, "geom"] <- 1 # gv1[9,"y"] <- 51 v2 <- terra::vect(gv1, "polygons") # plot(v2) # v2 <- is.valid(v2) terra::crs(v2) <- terra::crs(v) # v2 <- terra::makeValid(v2) if (getRversion() < "4.3.0") { # this same error crashes the session in R 4.3.0 when it is R-devel t10 <- try(postProcessTo(xVect, v2)) ## Error : TopologyException: Input geom 1 is invalid: ## Self-intersection at 6.0905735768254896 49.981782482072084 expect_true(!is(t10, "try-error")) } # Projection --> BAD BUG HERE ... CAN"T REPRODUCE ALWAYS --> use sf for testing Dec 9, 2022 if (FALSE) { utm <- terra::crs("epsg:23028") # sf::st_crs("epsg:23028")$wkt # albers <- sf::st_crs("epsg:5070")$wkt vutm <- terra::project(v, utm) } if (.requireNamespace("sf")) { # utm <- sf::st_crs("epsg:23028")#$wkt utm <- terra::crs("epsg:23028") # $wkt vsfutm <- sf::st_transform(vsf, utm) vutm <- terra::vect(vsfutm) res100 <- 100 rutm <- terra::rast(vutm, resolution = res100) rutm <- terra::rasterize(vutm, rutm, field = "NAME_2") vsfInUTMviaCRS <- postProcessTo(vsf, sf::st_crs(rutm)) expect_true(is(vsfInUTMviaCRS, "sf")) expect_true(terra::same.crs(vsfInUTMviaCRS, rutm)) # from is sf, to is SpatRast --> skip maskTo if (getRversion() >= "4.1" && isWindows()) { vsfInUTMviaSpatRast <- suppressWarningsSpecific( falseWarnings = "attribute variables are assumed", postProcessTo(vOrigsf, rutm) ) expect_true(is(vsfInUTMviaSpatRast, "sf")) expect_true(sf::st_crs(vsfInUTMviaSpatRast) == sf::st_crs(rutm)) expect_true(isTRUE(all.equal( round(terra::ext(rutm), 6), round(terra::ext(vsfInUTMviaSpatRast), 6) ))) } # Check for cases where `to` does not overlap with `from` ext1 <- terra::ext(vsf) + c(-2, 2, -2, 2) ext1SR <- terra::rast(ext1) if (.requireNamespace("raster")) { ext1Ra <- raster::raster(ext1SR) ext1Ex <- raster::extent(ext1Ra) } ext2 <- terra::vect(ext1) terra::crs(ext2) <- terra::crs(vsf) ext3 <- sf::st_as_sf(ext2) if (.requireNamespace("sf")) { expect_warning(expect_message(postProcessTo(vOrigsf, ext3))) # sf gives warning too expect_warning(expect_message(postProcessTo(terra::vect(vOrigsf), ext2)), "no intersection") # expect_warning(expect_error(postProcessTo(vOrigsf, ext3))) # sf gives warning too # expect_error(postProcessTo(terra::vect(vOrigsf), ext2)) } if (.requireNamespace("sf")) { if (.requireNamespace("sp")) { expect_error(postProcessTo(as(vOrigsf, "Spatial"), as(ext2, "Spatial"))) } if (.requireNamespace("raster")) { expect_error(postProcessTo(as(vOrigsf, "Spatial"), ext1Ra)) } } # if (Sys.info()["user"] %in% "emcintir") { # env <- new.env(parent = emptyenv()) # suppressWarnings( # b <- lapply(ls(), function(xx) if (isSpat(get(xx))) try(assign(xx, envir = env, terra::wrap(get(xx))))) # ) # save(list = ls(envir = env), envir = env, file = "~/tmp2.rda") # # load(file = "~/tmp2.rda") # # env <- environment() # # b <- lapply(ls(), function(xx) if (is(get(xx, env), "PackedSpatRaster") || is(get(xx, env), "PackedSpatVector")) try(assign(xx, envir = env, terra::unwrap(get(xx))))) # } t11 <- postProcessTo(elevRas, vutm) expect_true(terra::same.crs(t11, vutm)) # use raster dataset -- take the projectTo resolution, i.e., res100 t13 <- postProcessTo(elevRas, rutm) expect_true(identical(terra::res(t13)[1], res100)) expect_true(terra::same.crs(t13, vutm)) # no projection t12 <- postProcessTo(elevRas, cropTo = vutm, maskTo = vutm) expect_false(terra::same.crs(t12, vutm)) # projection with errors if (getRversion() >= "4.1" && isWindows()) { # bug in older `terra` that is not going to be fixed here utm <- terra::crs("epsg:23028") # This is same as above, but terra way if (getRversion() < "4.3.0") { # this same error crashes the session in R 4.3.0 when it is R-devel vutmErrors <- terra::project(v2, utm) mess <- capture_messages({ t13a <- postProcessTo(xVect, vutmErrors) }) ## Error : TopologyException: Input geom 1 is invalid: ## Self-intersection at 6095858.7074040668 6626138.068126983 expect_true(sum(grepl("error", mess)) %in% 1:2) # not sure why crop does not throw error in R >= 4.2 expect_true(sum(grepl("fixed", mess)) %in% 1:2) # not sure why crop does not throw error in R >= 4.2 expect_true(is(t13a, "SpatVector")) } else { v2 <- terra::makeValid(v2) vutmErrors <- terra::project(v2, utm) } # Switch from qs to rds with Cache opts <- options(reproducible.cacheSaveFormat = "qs") t13a <- Cache(postProcessTo(xVect, vutmErrors)) opts <- options(reproducible.cacheSaveFormat = "rds") t13a <- Cache(postProcessTo(xVect, vutmErrors)) opts <- options(reproducible.cacheSaveFormat = "qs") t13a <- try(Cache(postProcessTo(xVect, vutmErrors)), silent = TRUE) a <- try(ncol(t13a), silent = TRUE) expect_false(is(a, "try-error")) } # try NA to *To # Vectors vutmSF <- sf::st_as_sf(vutm) xVect2SF <- sf::st_as_sf(xVect2) t14 <- postProcessTo(xVect2, vutm, projectTo = NA) expect_true(terra::same.crs(t14, xVect2)) expect_false(terra::same.crs(t14, vutm)) if (getRversion() >= "4.1" && isWindows()) { # bug in older `terra` that is not going to be fixed here suppressWarningsSpecific( falseWarnings = "attribute variables", t14SF <- postProcessTo(xVect2SF, vutmSF, projectTo = NA) ) expect_true(terra::same.crs(t14SF, xVect2SF)) expect_false(terra::same.crs(t14SF, vutmSF)) } t15 <- postProcessTo(xVect2, vutm, maskTo = NA) expect_false(terra::same.crs(t15, xVect2)) expect_true(terra::same.crs(t15, vutm)) if (getRversion() >= "4.1" && isWindows()) { # bug in older `terra` that is not going to be fixed here suppressWarningsSpecific( falseWarnings = "attribute variables", t15SF <- postProcessTo(xVect2SF, vutmSF, maskTo = NA) ) expect_false(terra::same.crs(t15SF, xVect2SF)) expect_true(terra::same.crs(t15SF, vutmSF)) } t18 <- postProcessTo(xVect2, vutm, cropTo = NA) expect_false(terra::same.crs(t18, xVect2)) expect_true(terra::same.crs(t18, vutm)) suppressWarningsSpecific( falseWarnings = "attribute variables", t18SF <- postProcessTo(xVect2SF, vutmSF, cropTo = NA) ) expect_false(terra::same.crs(t18SF, xVect2SF)) expect_true(terra::same.crs(t18SF, vutmSF)) # Rasters t16 <- postProcessTo(elevRas, rutm, cropTo = NA) expect_false(terra::same.crs(t16, elevRas)) expect_true(terra::same.crs(t16, rutm)) expect_true(terra::ext(t16) >= terra::ext(rutm)) t17 <- postProcessTo(elevRas, rutm, projectTo = NA) expect_true(terra::same.crs(t17, elevRas)) expect_false(terra::same.crs(t17, rutm)) t19 <- postProcessTo(elevRas, rutm, maskTo = NA) expect_false(terra::same.crs(t19, elevRas)) expect_true(terra::same.crs(t19, vutm)) expect_true(sum(values2(t19), na.rm = TRUE) > sum(values2(t13), na.rm = TRUE)) # Raster with Vector t16 <- postProcessTo(elevRas, vutm, cropTo = NA) expect_false(terra::same.crs(t16, elevRas)) expect_true(terra::same.crs(t16, vutm)) t16SF <- postProcessTo(elevRas, vutmSF, cropTo = NA) expect_false(terra::same.crs(t16SF, elevRas)) expect_true(terra::same.crs(t16SF, vutmSF)) t17 <- postProcessTo(elevRas, vutm, projectTo = NA) expect_true(terra::same.crs(t17, elevRas)) expect_false(terra::same.crs(t17, vutm)) t17SF <- postProcessTo(elevRas, vutmSF, projectTo = NA) expect_true(terra::same.crs(t17SF, elevRas)) expect_false(terra::same.crs(t17SF, vutmSF)) t19 <- postProcessTo(elevRas, vutm, maskTo = NA) expect_false(terra::same.crs(t19, elevRas)) expect_true(terra::same.crs(t19, vutm)) # expect_true(sum(values2(t19), na.rm = TRUE) > sum(values2(t13), na.rm = TRUE)) t19SF <- postProcessTo(elevRas, vutmSF, maskTo = NA) expect_false(terra::same.crs(t19SF, elevRas)) expect_true(terra::same.crs(t19SF, vutmSF)) # expect_true(sum(values2(t19SF), na.rm = TRUE) > sum(values2(t13), na.rm = TRUE)) t21 <- postProcessTo(elevRas, projectTo = vutm) t21SF <- postProcessTo(elevRas, projectTo = vutmSF) t20 <- postProcessTo(elevRas, projectTo = terra::crs(vutm)) expect_true(all.equal(t20, t21)) expect_true(all.equal(t20, t21SF)) # This is now (Apr 24, 2023) because passing all to terra::project --> which keep square pixels expect_false(identical(terra::size(elevRas), terra::size(t20))) t20res250 <- postProcessTo(elevRas, projectTo = terra::crs(vutm), res = 250) expect_true(all(terra::res(t20res250) == 250)) ## same projection change resolution only (will likely affect extent) y2 <- terra::rast(crs = terra::crs(y), resolution = 0.008333333 * 2, extent = terra::ext(y)) y2 <- terra::setValues(y2, rep(1, terra::ncell(y2))) t22 <- postProcessTo(elevRas, to = y2, overwrite = TRUE) # not sure why need this; R devel on Winbuilder Nov 26, 2022 expect_true(terra::same.crs(t22, elevRas)) expect_true(terra::ext(t22) == terra::ext(y2)) ## "identical" may say FALSE (decimal plates?) expect_true(identical(terra::res(t22), terra::res(y2))) expect_false(identical(terra::res(t22), terra::res(elevRas))) vutmSF <- sf::st_as_sf(vutm) xVectSF <- sf::st_as_sf(xVect) ## It is a real warning about geometry stuff, but not relevant here err <- capture_error( # there is an error in R 4.0 Windows warn <- capture_warnings({ t22 <- postProcessTo(xVectSF, vutmSF) }) ) # } if (.requireNamespace("raster")) { ras1 <- raster::raster(tf5) ras2 <- raster::raster(tf6) ras1[] <- ras1[] ras2[] <- ras2[] r3 <- list(ras1, ras2) # this is for r3[[1]][] <- r3[[1]][] r3[[2]][] <- r3[[2]][] # bring to RAM # Raster & SpatVect ras1Small <- cropTo(ras1, t18) ras1SmallMasked <- maskTo(ras1, t18) # The warning is about some datum suppressWarnings(ras1SmallAll <- postProcessTo(ras1, t18)) expect_true(is(ras1Small, "Raster")) expect_true(is(ras1SmallMasked, "Raster")) expect_true(is(ras1SmallAll, "Raster")) # Check that extents are similar expect_true(terra::ext(ras1SmallAll) < terra::extend(terra::ext(t18), terra::res(ras1SmallAll) * 2)) # SpatRaster & Raster t20CroppedByRas <- cropTo(t20, ras1SmallAll) t20MaskedByRas <- maskTo(t20, ras1SmallAll) t20ProjectedByRas <- projectTo(t20, ras1SmallAll) t20AllByRas <- postProcessTo(t20, ras1SmallAll) # only does terra::rast 1x expect_true(is(t20AllByRas, "SpatRaster")) expect_true(is(t20CroppedByRas, "SpatRaster")) expect_true(is(t20MaskedByRas, "SpatRaster")) expect_true(is(t20ProjectedByRas, "SpatRaster")) expect_equal(terra::res(t20ProjectedByRas), terra::res(ras1SmallAll)) expect_equal(terra::res(t20AllByRas), terra::res(ras1SmallAll)) expect_equal(terra::ext(t20AllByRas), terra::ext(ras1SmallAll)) expect_equal(terra::ext(t20ProjectedByRas), terra::ext(ras1SmallAll)) expect_true( sum(!is.na(values2(ras1SmallAll))) < sum(!is.na(values2(t20ProjectedByRas))) ) expect_true( sum(!is.na(values2(ras1SmallAll))) == sum(!is.na(values2(t20AllByRas))) ) # The below was slightly off because RasterLayer was not exactly same as t20 proj spatRas1SmallAll <- projectTo(terra::rast(ras1SmallAll), t20) expect_true( # these are off b/c of projection probably abs(sum(!is.na(values2(spatRas1SmallAll))) - sum(!is.na(values2(t20MaskedByRas)))) <= 0 ) if (interactive()) { terra::plot(ras1SmallAll) terra::plot(t18, add = TRUE) terra::plot(t20AllByRas) terra::plot(t20CroppedByRas) terra::plot(t20MaskedByRas) terra::plot(t20ProjectedByRas) } w <- terra::vect("POLYGON ((0 -5, 10 0, 10 -10, 0 -5))") w1 <- terra::vect("POLYGON ((0 -5, 10 0, 10 -10, 4 -2, 0 -5))") w1a <- terra::vect("POLYGON ((20 15, 30 20, 30 10, 24 18, 20 15))") w2 <- rbind(w, w1, w, w1a, w) w3 <- fixErrorsIn(w2) expect_true(!all(terra::is.valid(w2))) expect_true(all(terra::is.valid(w3))) } } })