Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

rast() - using the win argument when reading in multiple files causes problems #964

Closed
dfriend21 opened this issue Jan 4, 2023 · 1 comment

Comments

@dfriend21
Copy link
Contributor

I've run into an issue where using the win argument when reading in several raster files into a single SpatRaster causes problems later on. See the example below.

library(terra)
#> terra 1.6.47

# save out 10 rasters, each to different files
ext <- ext(0, 10, 0, 10)
paths <- sapply(1:10, function(i){
  path_i <- tempfile(fileext = ".tif") # also tried .nc, .grd, and .asc - results were the same
  mat <- matrix(i, nrow = 10, ncol = 10)
  writeRaster(rast(mat, extent = ext), path_i)
  return(path_i)
})

# read them in, making use of the 'win' argument to specify a spatial subset
sub_ext <- ext(2, 6, 2, 6)
rasts <- rast(paths, win = sub_ext)

# check the results
rasts[[1]] # correct
#> class       : SpatRaster 
#> dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> window      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : file8ace2e40c17d.tif 
#> name        : lyr.1 
#> min value   :    >1 
#> max value   :    1<

rasts[[2]] # incorrect - has extent 'sub_ext' but dimensions 10 x 10 (and therefore resolution 0.4 x 0.4), which is incorrect
#> class       : SpatRaster 
#> dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
#> resolution  : 0.4, 0.4  (x, y)
#> window      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : file8ace1d5f0d44.tif 
#> name        : lyr.1 
#> min value   :    >2 
#> max value   :    2<

# try performing some operations on the rasters
rasts + 1          # works
#> class       : SpatRaster 
#> dimensions  : 4, 4, 10  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> names       : lyr.1, lyr.1, lyr.1, lyr.1, lyr.1, lyr.1, ... 
#> min values  :     2,     3,     4,     5,     6,     7, ... 
#> max values  :     2,     3,     4,     5,     6,     7, ...

rasts[[1:10]] + 1  # works
#> class       : SpatRaster 
#> dimensions  : 4, 4, 10  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> names       : lyr.1, lyr.1, lyr.1, lyr.1, lyr.1, lyr.1, ... 
#> min values  :     2,     3,     4,     5,     6,     7, ... 
#> max values  :     2,     3,     4,     5,     6,     7, ...

rasts[[1]] + 1     # works
#> class       : SpatRaster 
#> dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> name        : lyr.1 
#> min value   :     2 
#> max value   :     2

rasts[[2]] + 1     # doesn't work
#> Warning: /private/var/folders/w5/21t55k2s2l5g7kj4vvvwhshm0000gn/T/Rtmpk9UmqM/file8ace1d5f0d44.tif:
#> Access window out of range in RasterIO(). Requested (2,4) of size 10x10 on
#> raster of 10x10. (GDAL error 5)
#> class       : SpatRaster 
#> dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
#> resolution  : 0.4, 0.4  (x, y)
#> extent      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84

(rasts + 1)[[2]]   # works
#> class       : SpatRaster 
#> dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> name        : lyr.1 
#> min value   :     3 
#> max value   :     3


#============================================

# now try the same thing, but save the rasters as a single, multi-layered .tif
rasts2 <- rast(lapply(1:10, function(i){
  mat <- matrix(i, nrow = 10, ncol = 10)
  return(rast(mat, extent = ext))
}))
path2 <- tempfile(fileext = ".tif")
writeRaster(rasts2, path2)

rasts2a <- rast(path2, win = sub_ext)
rasts2a[[1]]      # correct
#> class       : SpatRaster 
#> dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> window      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : file8ace26cf9f01.tif 
#> name        : lyr.1 
#> min value   :    >1 
#> max value   :    1<

rasts2a[[2]]      # correct
#> class       : SpatRaster 
#> dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> window      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : file8ace26cf9f01.tif 
#> name        : lyr.1 
#> min value   :    >2 
#> max value   :    2<

rasts2a[[2]] + 1  # correct
#> class       : SpatRaster 
#> dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 2, 6, 2, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> name        : lyr.1 
#> min value   :     3 
#> max value   :     3

# this time everything works

Created on 2023-01-04 with reprex v2.0.2

Session info
─ Session info ───────────────────────────────────────────────────────
 setting  value
 version  R version 4.2.2 (2022-10-31)
 os       macOS Monterey 12.4
 system   aarch64, darwin20
 ui       RStudio
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Chicago
 date     2023-01-04
 rstudio  2022.12.0+353 Elsbeth Geranium (desktop)
 pandoc   NA

─ Packages ───────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 bench         1.1.2   2021-11-30 [1] CRAN (R 4.2.0)
 cli           3.4.1   2022-09-23 [1] CRAN (R 4.2.0)
 codetools     0.2-18  2020-11-04 [1] CRAN (R 4.2.0)
 Rcpp          1.0.9   2022-07-08 [1] CRAN (R 4.2.0)
 rlang         1.0.6   2022-09-24 [1] CRAN (R 4.2.0)
 rstudioapi    0.14    2022-08-22 [1] CRAN (R 4.2.0)
 sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.2.0)
 terra       * 1.6-47  2022-12-02 [1] CRAN (R 4.2.0)

 [1] /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library

──────────────────────────────────────────────────────────────────────
@rhijmans
Copy link
Member

rhijmans commented Jan 5, 2023

Thank you very much. I think this has now been fixed. I get

dim(x[[1]])
#[1] 4 4 1
dim(x[[2]])
#[1] 4 4 1
res(x[[1]])
#[1] 1 1
res(x[[2]])
[1] 1 1

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants