inlabru-org / inlabru

inlabru
https://inlabru-org.github.io/inlabru/
78 stars 21 forks source link

Error with no CRS and ipoints #89

Closed ASeatonSpatial closed 3 years ago

ASeatonSpatial commented 3 years ago

I couldn't create a reproducible example, something seems strange with this specific situation and I can't pin down what.

This works using the gorillas data removing crs information from the mesh and the sp polygons object and then creating ipoints:

library(inlabru)
library(raster)   # has function "crs<-" 
data(gorillas)

mesh = gorillas$mesh
g_poly = gorillas$boundary

crs(g_poly) = NA
mesh$crs = NULL 

g_poly_ip = ipoints(g_poly, mesh)    # runs no problem

On my example I get the following error

> sp::proj4string(study_area)
[1] NA
> mesh$crs
NULL
> ips = ipoints(region = study_area, domain = mesh)
Error in slot(x, "proj4string") : 
  no slot of name "proj4string" for this object of class "CRS"

Things get weirder... if I do raster::crs(study_area) before I run ipoints() then it works. I couldn't quite believe this, I've done it 4 or 5 times both ways and it never works without it and always works with it.

> raster::crs(study_area)
CRS arguments: NA 
> mesh$crs
NULL
> ips = ipoints(region = study_area, domain = mesh)
> 

I really don't get this, crs(study_area) should just be retrieving the crs information, not changing anything.

Traceback when it doesn't work is

> traceback()
6: slot(x, "proj4string")
5: comment(slot(x, "proj4string"))
4: is.na.CRS(crs)
3: is.na(crs)
2: fm_crs_is_null(domain_crs)
1: ipoints(region = study_area, domain = mesh)

Running debug I think the issue is domain_crs <- fm_ensure_crs(domain$crs).

From playing around in debug:

Browse[2]> fm_crs_is_null(domain$crs)
[1] TRUE
Browse[2]> fm_crs_is_null(domain_crs)
Error in slot(x, "proj4string") : 
  no slot of name "proj4string" for this object of class "CRS"

It seems fm_ensure_crs(domain$crs) does not return a CRS object that fm_crs_is_null() likes.

Some more info on domain_crs:

Browse[2]> domain_crs
CRS arguments: NA 
Browse[2]> str(domain_crs)
Formal class 'CRS' [package "sp"] with 1 slot
  ..@ projargs: chr NA
Browse[2]> class(domain_crs)
[1] "CRS"
attr(,"package")
[1] "sp"

I'm stumped from here. This doesn't happen when I debug on the gorillas data and I can't figure out what is funny about my example here. I think this issue will happen whenever fm_ensure_crs() is called since it returns sp::CRS(NA_character_) and this throws an error:

Browse[2]> fm_crs_is_null(sp::CRS(NA_character_))
Error in slot(x, "proj4string") : 
  no slot of name "proj4string" for this object of class "CRS"

The only connection I can see with the weird raster::crs() thing is that the output CRS arguments: NA appears both there and in domain_crs. I have the mesh and polygon saved as .RDS files, they are not big, I can provide them if that would be helpful.

> sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

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

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

other attached packages:
[1] rgeos_0.5-5        cowplot_1.1.0      scales_1.1.1       inlabru_2.1.14.900 ggplot2_3.3.2     
[6] INLA_20.10.12-1    foreach_1.5.1      sp_1.4-4           Matrix_1.2-18     

loaded via a namespace (and not attached):
 [1] rstudioapi_0.11  magrittr_1.5     splines_4.0.3    tidyselect_1.1.0 munsell_0.5.0   
 [6] colorspace_1.4-1 lattice_0.20-41  R6_2.5.0         rlang_0.4.8      dplyr_1.0.2     
[11] tools_4.0.3      rgdal_1.5-18     grid_4.0.3       gtable_0.3.0     withr_2.3.0     
[16] iterators_1.0.13 ellipsis_0.3.1   yaml_2.2.1       tibble_3.0.4     lifecycle_0.2.0 
[21] crayon_1.3.4     purrr_0.3.4      vctrs_0.3.4      codetools_0.2-16 glue_1.4.2      
[26] compiler_4.0.3   pillar_1.4.6     generics_0.0.2   pkgconfig_2.0.3 
finnlindgren commented 3 years ago

The root cause of the is.na is a bug in the sp 1.4-4 package, for sp:::is.na.CRS, that I've now reported upstream and made a workaround in inlabru, commit dbcf8c2. However, with your rgdal and sp versions I wouldn't have expected inlabru to call is.na at all, and the other unpredictabilities are also a bit strange. Is it possible you have PROJ4? (PROJ4 is meant to still work, but I don't have access to a PROJ4 system to test things on). Check before and after the raster::crs() call what inlabru:::fm_has_PROJ6() returns (in a fresh R session).

finnlindgren commented 3 years ago

Confirmed bug in sp, and the inlabru workaround works, and the code was triggered by pre-PROJ6 on Ubuntu 18.04.