rspatial / terra

R package for spatial data handling https://rspatial.github.io/terra/reference/terra-package.html
GNU General Public License v3.0
533 stars 86 forks source link

Incorrect behavior in `resample` when raster has a window set #1433

Open dfriend21 opened 6 months ago

dfriend21 commented 6 months ago

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

kadyb commented 6 months ago

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
sandra-ab commented 5 months ago

I have noticed the same issue within a loop.

My "fix" is to remove terra temporary files after each iteration: tmpFiles(remove=TRUE)

owenssam1 commented 1 month ago

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.