incorrect handling of "flipped" rasters?

[itati01]

Hi,

I used a SpatRaster object for a zonal statistics with exact_extract() and terra::extract(exact=T) and found suspicious values in the first output (unexpected NA, values too low). The SpatRaster was created from a Erdas Imagine file which was created from a NetCDF.

When I create a RasterLayer, raster::raster() gives the warning "In .rasterFromGDAL(x, band = band, objecttype, ...): data seems flipped. Consider using: flip(x, direction='y')". After explicitly flipping the RasterLayer (and assigning the correct CRS), exact_extract returns plausible results, albeit not identical to terra::extract.

Note: As expected, the SpatRaster is plotted correctly but not the original RasterLayer.

sessionInfo()
R version 4.1.1 (2021-08-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252
[4] LC_NUMERIC=C                    LC_TIME=German_Germany.1252    

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

other attached packages:
[1] exactextractr_0.7.2 rgdal_1.5-27        sp_1.4-5            terra_1.4-11       
[5] dplyr_1.0.7         sf_1.0-3  

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.7          pillar_1.6.4        compiler_4.1.1      class_7.3-19       
 [5] tools_4.1.1         lifecycle_1.0.1     tibble_3.1.5        lattice_0.20-44    
 [9] pkgconfig_2.0.3     rlang_0.4.11        DBI_1.1.1           e1071_1.7-9        
[13] exactextractr_0.7.2 s2_1.0.7            raster_3.5-2        generics_0.1.1     
[17] vctrs_0.3.8         hms_1.1.1           classInt_0.4-3      grid_4.1.1         
[21] tidyselect_1.1.1    glue_1.4.2          R6_2.5.1            fansi_0.5.0        
[25] tzdb_0.1.2          readr_2.0.2         purrr_0.3.4         magrittr_2.0.1     
[29] scales_1.1.1        codetools_0.2-18    ellipsis_0.3.2      units_0.7-2        
[33] assertthat_0.2.1    colorspace_2.0-2    utf8_1.2.2          KernSmooth_2.23-20 
[37] proxy_0.4-26        wk_0.5.0            munsell_0.5.0       crayon_1.4.1  

[dbaston]

exact_extract is using methods from raster or terra to read values from the underlying rasters. If raster is requiring you to call flip, then using exact_extract isn't going to exempt you from that. You can also call exact_extract with a terra SpatRaster directly, and I would recommend doing so for performance reasons alone.

[itati01]

Thanks for the reply. Without raster's warning I would not have been aware of the issue because terra was not giving any hint on flipping.

For me, the problem with SpatRaster is that terra::extract and plot are working as expected (cf. https://github.com/rspatial/terra/issues/25#issue-593742762), unfortunately unlike exact_extract. So, I assumed that something might be wrong here, or missing.

[dbaston]

Sorry, I'm not understanding the issue here. I would need code and/or data showing what is not working how you expect.

[itati01]

Sure. I provide some adjusted code from my script as well as a test grid and two test polygons in Test.zip.

# value grid
g <- "rr_annual_1995-2021.2000.img"
# zones as polygons
p <- read_sf("test.gpkg")
# no warning
grd.terra <- terra::rast(g)
# SpatRaster not flipped
plot(grd.terra)
# warning: flipped raster
grd.raster <- raster::raster(g) 
# is indeed flipped
plot(grd.raster)
# the plot of the flipped RasterLayer looks fine
grd.flipped <- raster::flip(grd.raster)
plot(grd.flipped)

# exact extract with terra
zstat.terra <- data.frame(terra::extract(grd.terra, terra::vect(p), fun="mean", exact=T, cells=T, layer="NUTS_ID", na.rm=T))
# exact_extract with SpatRaster: this is my issue
zstat.wrong <- exactextractr::exact_extract(grd.terra, p, fun="mean", force_df=T, append_cols="NUTS_ID")
# exact_extract with the original RasterLayer
zstat.wrong1 <- exactextractr::exact_extract(grd.raster, p, fun="mean", force_df=T, append_cols="NUTS_ID")
# exact_extract with the flipped RasterLayer
zstat.flip <- exactextractr::exact_extract(grd.flipped, p, fun="mean", force_df=T, append_cols="NUTS_ID")

zstat.terra
#  ID       rr
# 1  1 947.6912
# 2  2 865.1431
zstat.wrong
#   NUTS_ID     mean
# 1   AL012 631.5314
# 2   AL013 660.4338
zstat.wrong1
#   NUTS_ID     mean
# 1   AL012 640.8984
# 2   AL013 655.8774
zstat.flip
#   NUTS_ID     mean
# 1   AL012 960.4647
# 2   AL013 838.2286

Converting grd.terra to a RasterLayer gives again the warning while the SpatRaster from the original RasterLayer is plotted correctly. Does this mean that terra is not actually modifing the raster but considers the transformation otherwise? This could explain that grd.terra is identical to grd.raster but not to grd.flipped and also explain my issue with exact_extract. Despite the identity of grd.terra and grd.raster, zstat.wrong and zstat.wrong1 are different (the difference between terra::extract(exact=T) and the last exact_extract could be explained by different approaches).

raster::raster(grd.terra)
plot(terra::rast(grd.raster))
# grd.terra and g
rd.raster are identical (same extent, difference is zero)
terra::summary(grd.terra - terra::rast(grd.raster))
# grd.terra and grd.flipped have different extents
grd.terra - terra::rast(grd.flipped)
# Error: [-] extents do not match

[dbaston]

This makes the issue very clear, thank you.

[dbaston]

Using the fix recently committed to terra, this example now gives a correct answer for zstat.wrong. It still gives a different answer for zstat.flip because, as far as I can tell, raster does not correctly read the extent of this file.