Unexpected behaviour when extracting from raster with simple feature collection

[luigidolcetti]

Hi,

I am recently (sure in the past it was working correctly) experiencing a very weird behavior of exact_extract when extracting pixel from a raster using a simple feature collection. Essentially, let's say I have a collection with 2 rows/polygons, doing:

exact_extract(one_raster, one_sf_collection) ##### returns wrong values.

doing:

exact_extract(one_raster, one_sf_collection[1,]) #### just one row but the class of the object is the same, returns (apparently) correct values.

The same happens applying an aggregative function (e.g. mean).

Thank you in advance for any suggestion,

Luigi

PS, Here is my session info:

R version 4.1.0 (2021-05-18) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale: [1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.1252

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

other attached packages: [1] RUNIMC001_1.0.0 shinyBS_0.61.1 shiny_1.7.2 ncdf4_1.19

loaded via a namespace (and not attached): [1] Rcpp_1.0.9 lattice_0.20-44 class_7.3-19 assertthat_0.2.1 digest_0.6.29
[6] utf8_1.2.2 mime_0.12 plyr_1.8.7 R6_2.5.1 e1071_1.7-11
[11] pillar_1.8.1 rlang_1.0.6 rstudioapi_0.14 fontawesome_0.3.0 raster_3.6-3
[16] jquerylib_0.1.4 DT_0.25 exactextractr_0.9.1 textshaping_0.3.6 shinyjs_2.1.0
[21] htmlwidgets_1.5.4 munsell_0.5.0 proxy_0.4-27 compiler_4.1.0 httpuv_1.6.6
[26] pkgconfig_2.0.3 systemfonts_1.0.4 htmltools_0.5.3 tidyselect_1.1.2 tibble_3.1.8
[31] codetools_0.2-18 randomForest_4.7-1.1 fansi_1.0.3 dplyr_1.0.10 withr_2.5.0
[36] later_1.3.0 sf_1.0-8 keys_0.1.1 grid_4.1.0 lwgeom_0.2-8
[41] jsonlite_1.8.0 xtable_1.8-4 lifecycle_1.0.2 DBI_1.1.3 magrittr_2.0.3
[46] units_0.8-0 scales_1.2.1 KernSmooth_2.23-20 cli_3.4.1 cachem_1.0.6
[51] promises_1.2.0.1 sp_1.5-0 bslib_0.4.0 ellipsis_0.3.2 ragg_1.2.2
[56] generics_0.1.3 vctrs_0.4.1 tools_4.1.0 glue_1.6.2 purrr_0.3.4
[61] crosstalk_1.2.0 yaml_2.3.5 fastmap_1.1.0 colorspace_2.0-3 terra_1.6-17
[66] shinydashboard_0.7.2 classInt_0.4-7 memoise_2.0.1 sass_0.4.2

[dbaston]

I'm not seeing this with the test data included in the package. For example:

rast_fname <- system.file(file.path('sao_miguel', 'clc2018_v2020_20u1.tif'),
                          package = 'exactextractr')

poly_fname <- system.file(file.path('sao_miguel', 'concelhos.gpkg'),
                          package = 'exactextractr')

r <- terra::rast(rast_fname)
polys <- st_read(poly_fname, quiet = TRUE)

exact_extract(r, polys[1,], 'mean')
# [1] 18.49499
exact_extract(r, polys, 'mean')[1]
# [1] 18.49499

Can you share some data that I can use to reproduce?

[luigidolcetti]

Hi Dan, thank you for the super prompt replay. Please, have a look .... test.zip

thnx

[luigidolcetti]

Update.... loading the raster with raster::raster() so that the origin is on disk gives the problem, loading the raster all in memory with raster::readAll() the problem is solved (unfortunatelly I cannot have the whole bunch of data I have to analyse in mem, the file in the test.zip is just a crop).... but I am sure this problem wasn't there one year ago.

[luigidolcetti]

Update..... ok, apparently re-saving the raster file as tif it works fine.... so the problem must be with ncdf format.

[dbaston]

If I load the raster with terra::rast, then the results are consistent whether I am summarizing one or both features.

If I load the raster with raster, then I indeed see a different result if I am summarizing both features. I also see this warning:

Warning in x@ptr$crop(y@ptr, snap[1], extend[1], opt) :
  GDAL Message 1: dimension #1 (easting) is not a Longitude/X dimension.
Warning in x@ptr$crop(y@ptr, snap[1], extend[1], opt) :
  GDAL Message 1: dimension #0 (northing) is not a Latitude/Y dimension.

For summarizing multiple features, exactextractr may switch to terra behind the scenes. I wonder if that is causing incorrect values in this case due to GDAL's confusion about the X and Y axes.

[luigidolcetti]

I see, that would make sense (also temporarely since terra was released 1 or 2 years ago), infact I noticed that loading with terra::rast() it flips the raster.... mmmm that's annoying, I suppose those warnings may come from the fact the crs is set to na. Anyway thanks, Luigi

[dbaston]

If you think GDAL should be able to read this file correctly, you might want to open an issue there. I don't know the details of the CF conventions off the top of my head.