terra
terra copied to clipboard
Incorrect behavior in `resample` when raster has a window set
I've come across some incorrect behavior in resample()
that seems to be occurring when the raster being resampled has a window set. Running resample()
the first time works and produces the expected result. However, rerunning the code produces a raster without values. Here's a reproducible example:
library(terra)
#> terra 1.7.73
Sys.info()[1:2]
#> sysname release
#> "Windows" "10 x64"
ext <- ext(6, 6.2, 49.6, 49.8)
elev <- rast(system.file("ex/elev.tif", package = "terra"))
elev_sm <- rast(sources(elev), win = ext)
# works
rs <- resample(elev_sm, elev)
hasValues(rs)
#> [1] TRUE
# doesn't work
rs <- resample(elev_sm, elev)
hasValues(rs)
#> [1] FALSE
Created on 2024-02-19 with reprex v2.0.2
BTW: This seems to work if elev_sm
is in memory.
ext <- ext(6, 6.2, 49.6, 49.8)
elev <- rast(system.file("ex/elev.tif", package = "terra"))
elev_sm <- rast(sources(elev), win = ext)
set.values(elev_sm)
rs <- resample(elev_sm, elev)
hasValues(rs)
#> [1] TRUE
rs <- resample(elev_sm, elev)
hasValues(rs)
#> [1] TRUE
I have noticed the same issue within a loop.
My "fix" is to remove terra temporary files after each iteration:
tmpFiles(remove=TRUE)
I am having exactly the same issue where I did not have it before. When I reach the resample step in a loop of edits to a raster, it erases the data values on the raster. If you use terra::hasValues()
right after the resample step, it is FALSE. I think this is a bug in a newer version of terra because I have not run this code for about a year now. I have no idea how to fix this. I tried the suggestion from @sandra-ab, but it seemed to erase some essential data for terra to work in my session.
Here is the code I used that would not work. It loops through a set of rasters to format them and usually makes it past the 1st raster, only to fail on the 2nd:
for(a in seq_along(env.files)) {
# load each raster into temp object
rast.hold <- terra::rast(env.files[a])
# begin edits
# crop new rasters to extent
rast.hold <- terra::crop(x = rast.hold, y = ext.obj, overwrite = FALSE)
# mask the layers
rast.hold <- terra::mask(x = rast.hold, mask = mask_layer)
# set crs
rast.hold <- terra::project(x = rast.hold, y = "EPSG:4326", method = "bilinear")
# set origin
terra::origin(rast.hold) <- c(-0.0001396088, -0.0001392488)
#resample to fit the extent/resolution of the main layer
#use bilinear interpolation, since values are continuous
rast.hold <- terra::resample(x = rast.hold, y = main_layer, method = "bilinear")
# crop again to ensure aggregation doesnt increase extent
rast.hold <- terra::crop(x = rast.hold, y = ext.obj, overwrite = FALSE)
#write out the new resampled rasters!
terra::writeRaster(
x = rast.hold,
filename = file.path(mypath, "v1", paste0(output.files[a], "_global", ".tif")),
filetype = "GTiff",
overwrite = FALSE
)
# remove object once its done
rm(rast.hold)
}
UPDATE: I reverted terra by 2 versions (I went to 1.7-65), and that fixed the issue. I had been using version 1.7-78.