Serra314 commented 1 year ago

Issue CRS

Francesco Serafini 2022-09-09

I have an issue setting up the correct CRS. I found a workaround but I was wondering if this may be relevant for future development. The issue presents when I create a polygon and use that polygon as boundary to create a mesh and then, use the created mesh to fit an LGCP model. I am using R 4.2.1, INLA 22.9.2 and Inlabru 2.5.3.9001.

The polygon (which represents California) is defined as follows:

library(rgdal)
library(INLA)
library(inlabru)

Poly.coords <- data.frame(x = c(-125.2, -119.0, -119.0, -114.0, -113.1, -113.5, -113.6, -114.5, -117.1,-117.9,-118.4,-121.0,-121.6,-123.8, -125.4, -125.4),
                          y = c(43.0, 43.0, 39.4, 35.7, 34.3, 32.9, 32.2, 31.7, 31.5, 31.9, 32.8, 33.7, 34.2, 37.7, 40.2, 40.5))

I set the CRS this way and then try to create a mesh

crs_obj <- CRS(SRS_string='EPSG:4326')
Poly.bdy <- SpatialPolygons((list( Polygons(list(Polygon(Poly.coords)), ID="A1"))),  
                            proj4string = crs_obj)
mesh <- inla.mesh.2d(boundary=Poly.bdy,
                      max.edge=c(0.5, 1),
                      crs=crs_obj)

This produces the error

Error in inla.wkt_get_ellipsoid_radius(inla.crs_get_wkt(crs)) : Ellipsoid settings not found

Changing the CRS to the following I am able to create the mesh

crs_obj <- CRS(SRS_string='EPSG:3309')
Poly.bdy <- SpatialPolygons((list( Polygons(list(Polygon(Poly.coords)), ID="A1"))),  
                            proj4string = crs_obj)

mesh <- inla.mesh.2d(boundary=Poly.bdy,
                      max.edge=c(0.5, 1),
                      crs=crs_obj)

but not to fit the model because running the function

ipoints(samplers = Poly.bdy, domain = mesh)

I got the error

Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘coordinates’ for signature ‘“NULL”’

The only way I found to set the CRS in a way that is working is

crs_obj <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

With which everything works.

I was wondering if there is an issue on how I create the CRS object (which I suspect) or if the problem is more foundamental on how CRS are handled internally.

Thank you

finnlindgren commented 1 year ago

This is a bug in how inla.wkt_get_ellipsoid_radius() and fm_wkt_get_ellipsoid_radius() locate the ellipsoid information. For epsg:4326 it's located in a different type of tree structure than for the other crs, and the code logic doesn't know how to find it. The code expects "DATUM" containing "ELLIPSOID" but instead the 4326 crs has "ENSEMBLE" containing "ELLIPSOID".

It looks like a relatively easy fix but it also affects fm_wkt_set_ellipsoid_radius() and the corresponding INLA methods.

finnlindgren commented 1 year ago

I've made the corresponding update in the INLA code as well, so it should work after the next binary INLA build.

Serra314 commented 1 year ago

That sounds great! Does this mean that it should be available very soon?

Cheers,

Francesco

Il ven 9 set 2022, 18:01 Finn Lindgren @.***> ha scritto:

I've made the corresponding update in the INLA code as well, so it should work after the next binary INLA build.

— Reply to this email directly, view it on GitHub https://github.com/inlabru-org/inlabru/issues/154#issuecomment-1242161862, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFYE6R5TBZ7C4EWHL3H7KWDV5NNNXANCNFSM6AAAAAAQIVCDOE . You are receiving this because you authored the thread.Message ID: @.***>

finnlindgren commented 1 year ago

It depends on when H makes a new build; if you need it soon, I suggest email Haavard to explain.

inlabru-org / inlabru

Problem with defining CRS correctly #154

Issue CRS