inlabru-org / inlabru

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

Support integration of spherical and globe meshes #166

Closed finnlindgren closed 1 year ago

finnlindgren commented 2 years ago

From https://groups.google.com/g/r-inla-discussion-group/c/1f6Lj6ZbbEc/m/AtpuVhGxAgAJ :

Our samplers is a SpatialPolygonsDataFrame (called “sp”) in lat-lon coordinates, and the domain is an inla.mesh object (called mesh1) on the sphere (S^2) in Cartesian geocentric xyz coordinates. The reason that our mesh is in Cartesian xyz is we are modelling global processes, so our mesh has to be on S^2. We find we have to supply loc= Cartesian, geodetic xyz coordinates to inla.mesh.2d() to make meshes with full global coverage (i.e. to have longitude-wrapping and equal area). The inla.mesh.2d function seems to project lon-lat locations onto R^2.

When I ran ipoints(samplers=sp, domain=mesh1, group=”index”), I got the error

sp::over(samplers, integ_sp, returnList = TRUE) : identicalCRS(x, y) is not TRUE

Then I tried to use inla.spTransform(mesh1, inla.CRS("longlat")) to make mesh1 have the same CRS as the SpatialPolygons, but I got the error 'crs0' is an invalid coordinate reference object.

Part of reply:

The ipoints function doesn't have full support for global meshes yet, as it relies on conversion to a flat equal-area projection to compute the integration weights. I hope to soon have a new ipoints version to properly support this use case.

It's possible this issue can't be fully solved until we get more of the more general sf support working.

A potential workaround would be to make an integration scheme for the entire mesh, with the mesh nodes as integration points; something like this:

ips_sphere <- SpatialPointsDataFrame(mesh$loc, data=data.frame(
    weight = inla.mesh.fem(mesh)$va * earth.radius^2
  ), crs = fm_CRS("sphere"))
ips_globe <- fm_spTransform(ips_sphere, fm_CRS("globe"))

and then filter out the points outside the polygon, with a direct call to sp::over. That way, the mesh wouldn't need to be involved in the coordinate transformations.

For full integration accuracy along the sampler boundary, this step needs to be done before aggregating the integration points, but it might be sufficient as a temporary workaround.

The full version of the "workaround" is also the likely full solution, so something might be doable in the short term as well.