rspatial / terra

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

Using the win argument when reading in multiple files results in: [rast] extents do not match #1589

Closed Rapsodia86 closed 4 days ago

Rapsodia86 commented 2 months ago

Hi Robert! I am not sure if that is an issue, or the expected behavior. If the latter, then I take it back! I do have many files that I want to load as a raster stack. They have slightly different extents. To handle that, I would use the "win" argument in rast(), to "internally crop" the data and load up.

If I do this in a loop, everything is ok. But if I use a vector of file names, then the error shows up: Error: [rast] extents do not match.

Does that mean rast() first stacks all the rasters and then applies the "win"? Or it should be that it takes each object separately (like in a loop), applies the "win" and then stacks?

Below is a small example with only two files: (terra dev version)

library(terra)
terra 1.7.79
> aoi_ext <- ext(c(633138, 634573, 4694581, 4697348))
> rasts <- c("X:/MYDIR/HLS_S30_sr_evi2_doy2023342_aid0001_16N.tif", "X:/MYDIR/HLS_S30_sr_evi2_doy2023347_aid0001_16N.tif")
> 
> rstack <- rast(rasts,win=aoi_ext)
Error: [rast] extents do not match

##Here one by one works!:

> rstack <- rast()
> for (ff in rasts){
+   rstack <- c(rstack,rast(ff,win=aoi_ext))
+ }
Warning message:
[rast] the first raster was empty and was ignored 
> rstack
class       : SpatRaster 
dimensions  : 92, 47, 2  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
window      : 633150, 634560, 4694580, 4697340  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 16N (EPSG:32616) 
sources     : HLS_S30_sr_evi2_doy2023342_aid0001_16N.tif  
              HLS_S30_sr_evi2_doy2023347_aid0001_16N.tif  
names       : HLS_S30_sr_evi2~342_aid0001_16N, HLS_S30_sr_evi2~347_aid0001_16N 
sessionInfo()
R version 4.4.0 (2024-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] terra_1.7-79

loaded via a namespace (and not attached):
[1] compiler_4.4.0   tools_4.4.0      Rcpp_1.0.12      codetools_0.2-20

Thanks!

win_terra_rspatial_test.zip

rhijmans commented 4 days ago

This does not work as indicated by the error message

library(terra)
ff <- c("HLS_S30_sr_evi2_doy2023347_aid0001_16N.tif", "HLS_S30_sr_evi2_doy2023342_aid0001_16N.tif")
e <- ext(c(633138, 634573, 4694581, 4697348))
rast(ff)
# Error: [rast] extents do not match
rast(ff, win=e)
# Error: [rast] extents do not match

This is a work-around if you know what you are doing

x <- rast(lapply(ff, \(f) rast(f, win=e)))

But things would go wrong if you now do

window(x) <- NULL
x * 1
# Error: [*] too few values for writing: 90288 < 180576
#In addition: Warning message:
#C:/Users/rhijm/Downloads/win_terra_rspatial_test/HLS_S30_sr_evi2_doy2023342_aid0001_16N.tif: Access window out of #range in RasterIO().  Requested (0,0) of size 304x297 on raster of 304x278. (GDAL error 5) 

So to be on the safe side, I would remove your loophole --- but I won't do that for now.