coercing inla.mesh to Spatial* object

[bsp...@ncsu.edu]

Hi,

Is it possible to coerce an inla.mesh to a Spatial* object (e.g., SpatialPolygonsDataFrame)? I like to retain the polygons to use in other packages (e.g., tmap, ggplot2, etc).

Thanks!

[Finn Lindgren]

There is no SpatialTriangulation class unfortunately, so it would need to be a SpatialPolygons object. I never wrote such a converter, but we have implemented ggplot2 geoms for inla.mesh objects in the inlabru package; see the tutorials/workshop materials on inlabru.org for examples. For basic plotting of just the mesh,

ggplot() + gg(mesh)

should do it, but it can also map functions defined on the mesh vertices, etc.

I can see the benefit of having an explicit object conversion though. I’ll put it on my list...

Finn

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

[Finn Lindgren]

Hi again,

I've made a quick version of a mesh-to-Spatial* conversion function,
inla.mesh2sp. See code below.
It outputs SpatialPolygonsDataFrame (for the triangles) and
SpatialPoints (for just the vertices).
Including a special object for just the domain polygon is harder, but
we have internal code in the
'excursions' package that I could covert to work with this new function.
Is this the sort of thing you had in mind? Can you give it a try and
send me comments, and then I'll add it to the package.

Finn

#' Convert inla.mesh to sp objects
#'
#' @param mesh An \code{\link{inla.mesh}} object
#' @return A list with \code{sp} objects for triangles and vertices:
#' \describe{
#' \item{triangles}{\code{SpatialPolygonsDataFrame} object with the
triangles in
#' the same order as in the original mesh, but each triangle looping through
#' the vertices in clockwise order (\code{sp} standard) instead of
#' counterclockwise order (\code{inla.mesh} standard). The \code{data.frame}
#' contains the vertex indices for each triangle, which is needed to link to
#' functions defined on the vertices of the triangulation.
#' \item{vertices}{\code{SpatialPoints} object with the vertex coordinates,
#' in the same order as in the original mesh.}
#' }
#' @export
inla.mesh2sp <- function(mesh) {
crs <- inla.CRS(inla.CRSargs(mesh$crs))
isgeocentric <- identical(inla.as.list.CRS(crs)[["proj"]], "geocent")
if (isgeocentric || (mesh$manifold == "S2")) {
stop(paste0(
"'sp' doesn't support storing polygons in geocentric coordinates.\n",
"Convert to a map projection with inla.spTransform() before
calling inla.mesh2sp()."))
}

triangles <- SpatialPolygonsDataFrame(
Sr = SpatialPolygons(lapply(
1:nrow(mesh$graph$tv),
function(x) {
tv <- mesh$graph$tv[x, , drop = TRUE]
Polygons(list(Polygon(mesh$loc[tv[c(1, 3, 2, 1)],
1:2,
drop = FALSE])),
ID = x)
}
),
proj4string = crs
),
data = as.data.frame(mesh$graph$tv[, c(1, 3, 2), drop = FALSE]),
match.ID = FALSE
)
vertices <- SpatialPoints(mesh$loc[, 1:2, drop = FALSE], proj4string = crs)

list(triangles = triangles, vertices = vertices)
}

[bsp...@ncsu.edu]

Thanks, Finn, this seems to be working as I hoped. No current suggestions for improvements but I will let you know as I continue to work with it. Thank you for a speedy response!

[Jenesais Pas]

Thank you very much. This function is really great. It also helps to export the mesh into a Desktop GIS (e.g. QGIS for mapping purpose). I am wondering if the function is in the meanwhile already implemented in any packages?

[Finn Lindgren]

Hi,

I had forgotten about this function! Thanks for the reminder...

I'll add it to the INLA package soon (Haavard is transferring the code
mercurial repository to git, and we just need to make sure that the
transfer worked properly before making further changes, and then
there's also an R3-R4 transition needed),
I'll also update it to work with newer CRS methods (anyone with a new
PROJ6 installation will have noticed that rgdal/sp now gives lots of
warnings for old CRS specifications). I'm working on a blog post about
this for INLA users.

Finn

> --
> You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/4338feb4-7fdc-49f5-8f37-04d7902518c9o%40googlegroups.com.



--
Finn Lindgren
email: finn.l...@gmail.com

[Jenesais Pas]

Hi,

thanks for the quick response and the update.

It would also be really great if the function is more generic, so that it can return a sp-object also for inla.nonconvex.hull. That enables the possibility to interactively look at intermediate results via, e.g., mapview. (Btw, I was wondering why all spatial objects are in sp-manner and not sf - as simple features are more flexible by, e.g., allowing geometrycollections of points and lines or polygons).

> To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsub...@googlegroups.com.

[Finn Lindgren]

Hi,

from what I recall, the reason I never wrote a complete
inla.segment-to-sp converter is the handling of holes; the
inla.segment objects doesn't explicitly label edges as holes or not;
it only keeps track of what is on the inside and outside of each edge
segment. To identify holes, it would need to build a mesh, which
_does_ correctly identify the emergent holes from the segment
representation. It's not in principle difficult, I just never got
around to it.

Similarly, sf did not yet exist when the sp support was added to INLA,
and there hasn't been time to add sf support since when sf appeared;
it would touch on so much of the mesh related code that it would be
better to add that support properly as part of converting the fmesher
interface to a separate package, which is on my TODO-list!

Finn

On Thu, 2 Jul 2020 at 09:51, 'Jenesais Pas' via R-inla discussion

>> > To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.

>> > To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/4338feb4-7fdc-49f5-8f37-04d7902518c9o%40googlegroups.com.
>>
>>
>>
>> --
>> Finn Lindgren
>> email: finn.l...@gmail.com
>

> --
> You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/cc922362-fb15-4bb0-926c-f71df281aa4fo%40googlegroups.com.

[Russell Dinnage]

I have been getting lots of warnings working with the new version of sf and then converting to sp Spatial objects to work with INLA. It is great to hear you are working on updating to the new CRS methods. I hope this can resolve this for me. Also very excited to see INLA moving to github! 

Anyway, just for your information as you move towards updating to R 4 and new CRS methods, I thought I'd post my problem, in case it helps:

When using inla.mesh.2d with crs = INLA::inla.CRS("sphere"), and using sp data that was in "EPSG:4326" projection converted from sf, I get more than 50 warnings of this variety:

1: In showSRID(SRS_string, format = "PROJ", multiline = "NO") :

  Discarded ellps unknown in CRS definition: +proj=geocent +R=1 +units=m +no_defs

2: In showSRID(SRS_string, format = "PROJ", multiline = "NO") :

  Discarded datum unknown in CRS definition

3: In showSRID(SRS_string, format = "PROJ", multiline = "NO") :

  Discarded ellps unknown in CRS definition: +proj=geocent +R=1 +units=m +no_defs

4: In showSRID(SRS_string, format = "PROJ", multiline = "NO") :

  Discarded datum unknown in CRS definition

When I plot the mesh in rgl, it looks okay, however, if I use the mesh in INLA::inla.mesh.project, with some coordinates also converted to sp from sf with "EPSG:4326" projection I get the same kind of warnings, and then when I plot the points on the globe, they don't look quite right especially near the poles.

I'm stuck at this point! I don't know whether to trust the results, given the warnings, and the plot not looking right. This code used to work before I upgraded my sf version (and presumably with it my PROJ6 version). But I guess I will wait for that blog post you mentioned which will hopefully shed light on this! 

Thanks for all your work on INLA!

Cheers,

Russell

[Finn Lindgren]

Hi,

The blog post has been held up by teaching and other work, but the short answer is that the current version of inla does have proj6 support. If the function INLA::inla.has_PROj6() exists and returns TRUE you have a new INLA and proj6.

Even with the new version you will still receive lots of warnings (but different warnings than before) that come from the rgdal package, but most of these can be ignored, since they only relate to the backwards compatibility behaviour of rgdal for the old proj4 strings, that are now mostly ignored, and the actual projection info stored as a “WKT” string. The new inLa version has some tools to help work with those, but the bulk of the warnings can be turned off in rgdal. (And as I understand it they will deactivate these compatibility warnings by default in a future version)

Finn

On 6 Jul 2020, at 03:35, Russell Dinnage <r.di...@gmail.com> wrote:



--

You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/122f3140-3e45-413c-971a-192460e58aa2o%40googlegroups.com.

[Sylvan Benaksas]

Hi Finn,

I wonder if there has been any update on this, I am currently trying to convert mesh to polygons using your code from the earlier email but I am getting many warnings

Warning messages:
1: In inla.not_for_PROJ6("inla.CRSargs") :
  'inla.CRSargs' should not be used with PROJ6 and rgdal v3
Call stack for developer debugging:
inla.mesh2sp(mesh)
inla.CRS(inla.CRSargs(mesh$crs))
inla.CRSargs(mesh$crs)
inla.not_for_PROJ6("inla.CRSargs")
2: In inla.not_for_PROJ6("inla.as.list.CRS") :
  'inla.as.list.CRS' should not be used with PROJ6 and rgdal v3
Call stack for developer debugging:
inla.mesh2sp(mesh)
inla.as.list.CRS(crs)
inla.not_for_PROJ6("inla.as.list.CRS")
3: In inla.not_for_PROJ6("inla.as.list.CRSargs") :
  'inla.as.list.CRSargs' should not be used with PROJ6 and rgdal v3
Call stack for developer debugging:
inla.mesh2sp(mesh)
inla.as.list.CRS(crs)
inla.as.list.CRSargs(inla.CRSargs(x))
inla.not_for_PROJ6("inla.as.list.CRSargs")
4: In inla.not_for_PROJ6("inla.CRSargs") :
  'inla.CRSargs' should not be used with PROJ6 and rgdal v3
Call stack for developer debugging:
inla.mesh2sp(mesh)
inla.as.list.CRS(crs)
inla.as.list.CRSargs(inla.CRSargs(x))
inla.CRSargs(x)
inla.not_for_PROJ6("inla.CRSargs")

Cheers,

Sylvan

[Finn Lindgren]

The old code likely needs to be updated, but we're in the process of
implementing sf support since rgdal etc are being retired soon. In the
meantime, the fm_* functions in inlabru may have more bug fixes than
the sp-functions in INLA; eventually all of these will be in the new
fmesher package instead to make it easier to maintain a coherent
codebase. We've started an sf-support branch for inlabru that has an
experimental mesh-to-sf converter.
Whether that's useful depends on what problem you're trying to solve.
Finn

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/756d97e5-25b1-4bf3-8931-d9ec80ec8725n%40googlegroups.com.

[Sylvan Benaksas]

Hi Finn,

Thanks for the update, it has actually worked for me, as you said before the warnings can be ignored,

Cheers,

Sylvan