create mesh using polygon with open water

[cp...@terpmail.umd.edu]

I am trying to make a 2d inla mesh that will contain only water. As you can see below, my coastline is open at the edges because water keeps flowing off the edges of the map. Is it ok to use this as my boundary for making a mesh? Or is there some way to "seal off" those open water edges?

#find coastline boundaries by extracting raster cells that border NA values (NA values are land while all other values are water)

RasterBounds <- raster::boundaries(raster, type = "inner", classes = F, directions = 4)

#convert boundary cells to polygon object

pol <- rasterToPolygons(RasterBounds, fun=function(x){x==1})

plot(pol)

#create 2d mesh using polygon as a boundary, with max.edge equal to twice the size of 1 of my raster pixels (1 pixel=100mX100m)

mesh2 = inla.mesh.2d(boundary = pol,

                     max.edge = 200)

plot(mesh2)

[Finn Lindgren]

Hi,

If you can extract the coastline as “linetrains”, e.g. as a SpatialLines object, then all you need to “seal off” the domain is to add another set of lines, or a whole polygon, that leads to the combined shape being closed; the mesh will be built on the domain that the input lines&polygons lie “inside” all of them. For lines, that is determined by counterclockwise ordering, I.e. “inside” is to your left as you traverse the lines. I’m not on my laptop at the moment, but IIRC the syntax is

inla.mesh.2d(boundary=list(list(mylines, myenclosingpolygon)), max.edge=..., cutoff=...)

You’ll most likely need cutoff greater than zero to avoid tiny triangles in the fjords.

That said, a mesh of this type will have very strong boundary effects in the fjords.

A more appealing alternative from a modelling perspective is to use the model developed by Haakon Bakka, which has been mentioned on this list previously; there, the mesh stretches over land as well, but with such a short correlation range there that the end result is information is forced to flow in the ocean only, but without the strong boundary effects.

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.

[Haakon Bakka]

Hi,

Francesco Serafini and me are creating a tutorial for coastline polygons and meshes, and one for the Barrier model. I put the preliminary versions online because of your question:

https://haakonbakka.bitbucket.io/btopic127.html

https://haakonbakka.bitbucket.io/btopic128.html

These should be able to answer most of your questions. 

When it comes to cutoff, I use 1/5 of the max.edge.length argument. This works out well (almost always). But look at the mesh to see that it does not look weird in places.

Kind regards,

Haakon Bakka

Finn

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@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.

--

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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.

[Carrie Perkins]

Thank you both so much! I think I finally got it to work using a mixture of both of your code, plus guidance from this post: https://gis.stackexchange.com/questions/187798/create-polygons-of-the-extents-of-a-given-raster-in-r

Here is what I did:

#extracted extent from my raster

e <- extent(raster1)

#turned extent into SpatialPolygons object to create the overall 'study space' polygon

p <- as(e, 'SpatialPolygons')

#got a polygon that surrounds cells in raster that are not NA to create the 'water' polygon

#first made a raster containing only the water cells (excluded land cells, which were NAs)

r <- raster1 > -Inf

#then converted water raster to polygons (you need to have package 'rgeos' installed for this to work)

pp <- rasterToPolygons(r, dissolve = T)

#plot study space polygon

plot(p, lwd=5, border='red')

#plot water polygon on top of it

plot(pp, lwd=3, border='blue')

#create mesh, with pp (the water polygon) on the left because it is inside the study space polygon, p

mesh1 <- inla.mesh.2d(boundary=list(list(pp, p)), max.edge=c(1000,2000), cutoff=200)

#plot mesh

plot(mesh1)

I do have a few additional questions: each pixel in my original raster was 100m by 100m and I have covariate data for each pixel that I want to use in an R inla bayesian spatial regression model. Does that mean that there should be a triangle for each pixel? If yes, should I change the max edge to 100? 

Also, what determines the amount of land the mesh stretches over in the Haakon Bakka model? Is the short correlation range dictated by the cutoff being 1/5 of the max edge?

Is this mesh only used to help the model decide where to predict the results spatially, or does the extent of the mesh need to match the data points perfectly?

Thank you again!!

Finn

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

To post to this group, send email to r-inla-discussion-group@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.

--
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-group+unsubscr...@googlegroups.com.

[Haakon Bakka]

Hi,

See inline.

On 30 April 2018 at 20:42, Carrie Perkins <cp...@terpmail.umd.edu> wrote:

Thank you both so much! I think I finally got it to work using a mixture of both of your code, plus guidance from this post: https://gis.stackexchange.com/questions/187798/create-polygons-of-the-extents-of-a-given-raster-in-r

I do have a few additional questions: each pixel in my original raster was 100m by 100m and I have covariate data for each pixel that I want to use in an R inla bayesian spatial regression model. Does that mean that there should be a triangle for each pixel? If yes, should I change the max edge to 100? 

No. The triangles are only for representing the spatial random effect. You can still observe pixels. Just pretend they are locations and at each location you have covariate values.

 

Also, what determines the amount of land the mesh stretches over in the Haakon Bakka model? Is the short correlation range dictated by the cutoff being 1/5 of the max edge?

The cutoff is related to the approximation of the coastline. How far you stretch the mesh over land is the same as if you didi a stationary model.  What you have seems reasonable.

 

Is this mesh only used to help the model decide where to predict the results spatially, or does the extent of the mesh need to match the data points perfectly?

It is only used for the spatial random effect. You can only predict or fit data in the "inner" part of the mesh. Where exactly you have datapoints (or prediction points) is mostly irrelevant (as long as it is not outside the inner mesh).

 

Haakon

[Carrie Perkins]

Thank you both very much! I understand now.