Error when using read_RAST() with raster map containing values beyond int16

[TillF]

Thanks for providing rGRASS.

When trying to read a raster map with large values, read_RAST() fails:

execGRASS(cmd = "r.mapcalc", expression="test_t=66000", flags="overwrite") #create raster with values larger than int32
execGRASS("r.info", flags = "r", map = "test_t", intern = TRUE)
tt = read_RAST(vname = "test_t") #try to import raster: fails with error "[rast] file does not exist"

tt = read_RAST(vname = "test_t", NODATA=70000) #works

possible culprit within readRAST:

NODATAi=4294967295 #NA-value assigned within read_RAST(), probably too high?
execGRASS("r.out.gdal", input = "test_t", output = tempfile(),   
          nodata = NODATAi, flags = "overwrite") #fails

fails with

WARNUNG: Mismatch between metadata nodata value and actual nodata value in
         exported raster: specified nodata value 4.29497e+09 gets converted
         to -2147483648 by selected GDAL datatype.
WARNUNG: GDAL-Datentyp: Int32, g�ltiger Wertebereich: -2147483648 -
         2147483647
FEHLER: Export der Rasterkarte abgebrochen.

Using a smaller value for NODATAi works:

execGRASS("r.out.gdal", input = "test_t", output = tempfile(),  #works (with warnings) 
          nodata = 70000, flags = "overwrite")

I suggest revising the determination of NODATAi in read_RAST() or guiding the user to provide appropriate values for the argument NODATA with some helpful error message.

tested on GRASS 8.2.1 (2023)

sessionInfo() R version 4.2.2 (2022-10-31 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale: [1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8 LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.utf8

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

other attached packages: [1] lumpR_3.1.72 sf_1.0-9 rgeos_0.6-1 RODBC_1.3-20 foreach_1.5.2
[6] xts_0.12.2 zoo_1.8-11 maptools_1.1-8 cluster_2.1.4 rgrass_0.3-9
[11] XML_3.99-0.13 terra_1.7-55 remotes_2.4.2 pryr_0.1.6 devtools_2.4.5
[16] usethis_2.1.6 panelaggregation_0.1.1 data.table_1.14.8 curl_5.0.0 gWidgets_0.0-54
[21] rpart_4.1.19 rgdal_1.6-7 raster_3.6-26 sp_1.5-1

loaded via a namespace (and not attached): [1] Rcpp_1.0.11 lattice_0.20-45 class_7.3-20 prettyunits_1.1.1 ps_1.7.2 utf8_1.2.2 digest_0.6.31
[8] mime_0.12 R6_2.5.1 e1071_1.7-12 pillar_1.8.1 rlang_1.0.6 rstudioapi_0.14 miniUI_0.1.1.1
[15] callr_3.7.3 urlchecker_1.0.1 stringr_1.5.0 foreign_0.8-83 htmlwidgets_1.6.0 proxy_0.4-27 shiny_1.7.4
[22] compiler_4.2.2 httpuv_1.6.9 pkgconfig_2.0.3 pkgbuild_1.4.0 htmltools_0.5.4 tidyselect_1.2.0 tibble_3.1.8
[29] codetools_0.2-18 fansi_1.0.3 dplyr_1.0.10 crayon_1.5.2 later_1.3.0 grid_4.2.2 xtable_1.8-4
[36] lifecycle_1.0.3 DBI_1.1.3 magrittr_2.0.3 units_0.8-1 KernSmooth_2.23-20 cli_3.4.1 stringi_1.7.8
[43] cachem_1.0.6 fs_1.5.2 promises_1.2.0.1 generics_0.1.3 ellipsis_0.3.2 vctrs_0.5.1 iterators_1.0.14
[50] tools_4.2.2 glue_1.6.2 purrr_1.0.1 processx_3.8.0 pkgload_1.3.2 fastmap_1.1.0 sessioninfo_1.2.2 [57] classInt_0.4-8 memoise_2.0.1 profvis_0.3.7

[veroandreo]

I can reproduce the error, too.

GRASS 8.2.1

sessionInfo()
R version 4.2.3 (2023-03-15)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: Fedora Linux 37 (Xfce)

Matrix products: default
BLAS/LAPACK: /usr/lib64/libflexiblas.so.3.3

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
[1] terra_1.7-55  rgrass_0.3-8  XML_3.99-0.14

loaded via a namespace (and not attached):
[1] compiler_4.2.3    tools_4.2.3       rstudioapi_0.15.0 Rcpp_1.0.11       codetools_0.2-19 

[rsbivand]

@veroandreo I'll take a look, no obvious reason. The problem may be code trying to guess a NODATA value if not given.

If you could respond to the review request for PR #81 that would help, as no fix for the current problem will be published before we shift to xml2 from EOL XML.

[rsbivand]

@TillF Which test location are you using? And which version of GDAL was GRASS build against, on which platform (Windows)? I can reproduce the problem with Fedora 39, GRASS 8.3.1, GDAL 3.8.1.

[veroandreo]

@veroandreo I'll take a look, no obvious reason. The problem may be code trying to guess a NODATA value if not given.

If you could respond to the review request for PR #81 that would help, as no fix for the current problem will be published before we shift to xml2 from EOL XML.

yes, I have it bookmarked to do asap

[TillF]

I found this error when working with my own location with the following dimensions: projection: 1 (UTM) zone: 36 datum: wgs84 ellipsoid: wgs84 north: 1769869.81527969 south: 854105.24287058 west: 446314.28520378 east: 1246996.97157839 nsres: 30.95263207 ewres: 30.95263207 rows: 29586 cols: 25868 cells: 765330648

I am running Windows10, see details in sessionInfo() in my first post. I don't know how to get the GDAL version. Please let me know, if this is needed.

[rsbivand]

@TillF How did you install GRASS, the stand-alone installer or OSGeo? Both bundle GDAL, so knowing how GRASS was installed should indicate the version of GDAL.

[TillF]

I don't quite remember, but seeing it resides in C:/Program Files/GRASS GIS 8.3 I assume it's the stand-alone version.

[rsbivand]

OK, I've located the problem (perhaps), which may be a change in r.out.gdal or GDAL, previously the choice of UInt32 happened automatically based on the values and NODATA, now it looks at the values and gets to Int32, so NODATA is too big. I've added explicit declaration of UInt32 and UInt16 where the NODATA guesses require that. The commit is to the xml2 branch. I'll try to build a Windows binary package for testing, later.

[rsbivand]

This is the built source package from the xml2 branch with the NODATA fix: rgrass_0.4-1.tar.gz

[rsbivand]

And the Windows binary package, built with Rtools43 for R 4.3: rgrass_0.4-1.zip

[TillF]

Thanks for looking into that so swiftly. I updated with the Windows-binary. With tt = read_RAST(vname = "test_t") I now get set NODATA manually to a feasible value As this is difficult to understand for uninformed users like me, I suggest to change this to stop(paste0("Set NODATA manually to a feasible value, try e.g.", lres$max+1))

However, calling the function with argument NODATA, now causes an error: tt = read_RAST(vname = "test_t", NODATA=70000) #error: Objekt 'typei' nicht gefunden

[rsbivand]

OK, I'll look at this again. I'm open to dropping NODATA completely, forcing the user to specify it manually in all cases. Most of the guesses are wrong, I'm afraid.

[rsbivand]

New source package (moving the assignment of typei earlier): rgrass_0.4-1.tar.gz

[rsbivand]

New Windows binary package: rgrass_0.4-1.zip

[TillF]

Thanks, the second update solved the problem for my test case, i.e. even without having to supply the value for NODATA manually. However, the behaviour now is somewhat misleading, as trying to load a nonexistent raster map yields the newly-introduced error message:

> tt = read_RAST(vname = "does_not_exist") #
Fehler in read_RAST(vname = "does_not_exist") : 
  set NODATA manually to a feasible value

I reckon the function should better indicate that the respective raster was not found.

[rsbivand]

Thanks, will check. However, I cannot reproduce this:

> tt = read_RAST(vname = "does_not_exist")
Error in execGRASS("r.info", flags = "r", map = vname[i], intern = TRUE,  : 
  The command:
r.info -r map=does_not_exist
produced an error (1) during execution:
ERROR: Raster map <does_not_exist> not found

as r.info is called first to collect information on the raster to be read at: https://github.com/rsbivand/rgrass/blob/6611c3d304d91c3c7c918e72696b4bf1c2ce2904/R/rast_link.R#L67-L68. I'm on Linux, will check on Windows.

[rsbivand]

I can reproduce the problem on Windows, but see:

> tt = read_RAST(vname = "does_not_exist")
Error in read_RAST(vname = "does_not_exist") : 
  set NODATA manually to a feasible value
In addition: Warning messages:
1: In system(syscmd, intern = intern, ignore.stderr = ignore.stderr,  :
  running command 'r.info.exe -r map=does_not_exist' had status 1
2: In ifelse(x == "NULL", as.numeric(NA), as.numeric(x)) :
  NAs introduced by coercion
3: In system(syscmd, intern = intern, ignore.stderr = ignore.stderr,  :
  running command 'r.info.exe -g map=does_not_exist' had status 1

In Windows, a different system call is made, giving a warning not an error. Do you also see the warnings? I'll try to add extra logic to test for existence, and to go through NODATA setting only if there are NODATA cells.

[rsbivand]

New source package: rgrass_0.4-1.tar.gz

[TillF]

Yes, the warnings are vissible, I just failed to check them. Anyway, adding some extra logic in this "severe" case for raising a true error instead of a mere warning strikes me as a good idea, especially as the warnings are not explicitly stating that the map is missing.

[rsbivand]

With: rgrass_0.4-1.zip

>  tt = read_RAST(vname = "does_not_exist")
Error in read_RAST(vname = "does_not_exist") : does_not_exist not found

[TillF]

Perfect, works for me. Issue can be closed as far as I am concerned. Thanks.

[rsbivand]

Thanks!