SPDE model on S^2 / Unit of parameter spatial scale kappa

[gg]

Making simulation and inference under the SPDE model on S^2 with spatial (Lon,Lat) coordinates in degrees, 

in what unit is the scale parameter kappa expressed? 

[Finn Lindgren]

Do you mean for a model on S^2 embedded in R^3, or a model in
(Lon,Lat) as a subset of R^2?
I'm guessing you mean the former, in which case the answer is that in
r-inla, meshes on S^2 are by default constructed with radius 1. When
the true radius of the globe is given by R, say, the correlation range
then goes as R/kappa.

If given coordinates with a radius different than 1, inla.mesh.create
should produce the appropriately scaled mesh with radius not equal to
1, but this has to my knowledge not been used when modelling by
anyone.

Finn

[gg]

Many thanks for your prompt answer.

A related issue: 

I am interesting in getting a vector of points values from the SPDE Matérn model on S^2 

where the sampling sites belong to a small area on the sphere (large enough for the curvature to matter).

Borrowing from http://www.r-inla.org/examples/case-studies/simpson2011 

I got something along the line of the below which seems to work fine for a  window  large enough 

but returns many NAs for a small window. Any work around?

Many thanks beforehand!

Glles  

###############"

require(INLA);require(fields)

set.seed(1)

n <- 100

sigma2.approx = 0.5^2

range.approx = 0.4

kappa = sqrt(8)/range.approx

tau = 1/(4*pi*kappa^2*sigma2.approx)^0.5

## here sampling window smaller than needed to excerbate the problem

s <- cbind(runif(n,0,5),runif(n,40,45))

lonlat3D=function(lon,lat){

      cbind(cos((lon/180)*pi)*cos((lat/180)*pi),

            sin((lon/180)*pi)*cos((lat/180)*pi),

            sin((lat/180)*pi))

    }

true.radius.of.earth = 6371

s3 <- lonlat3D(s[,1],s[,2])

mesh <- inla.mesh.create(loc=s3,

                         cutoff=100/true.radius.of.earth,

                         refine=list(max.edge=500/true.radius.of.earth))

spde=inla.spde.create(mesh, model="matern")

## proj = inla.mesh.projector(mesh,dims=c(180,90),projection="longsinlat")

proj = inla.mesh.projector(mesh,dims=c(180,90),projection="longlat") ## better than longsinlat in terms of not getting NAs 

x <- inla.spde.query(spde, sample=list(tau=tau, kappa2=kappa^2))$sample

omesh = old.mesh.class(mesh)

cp=colorRampPalette(c("blue","cyan","yellow","red"))

plot(omesh, x, color.palette=cp, draw.edges=TRUE, draw.vertices=TRUE)

plotdata.x = inla.mesh.project(proj, x)

mean(is.na(plotdata.x))

dev.new();image.plot(plotdata.x, col=cp(256))

## extracting values from matrices for sampling sites

xdat <- numeric(nrow(s))

xmat <- matrix(nr=180,nc=90,plotdata.x)

for(k in 1:nrow(s))

  {

    xc <- s[k,1]

    yc <- s[k,2]

    ii <- which.min(abs(xc-seq(-180,180,length=180)))

    jj <- which.min(abs(yc-seq(-90,90,length=90)))

    xdat[k] <- xmat[ii,jj]

  }

mean(is.na(xdat)) ## many NAs

 

[Finn Lindgren]

On 19/03/13 14:43, gg wrote:
> Many thanks for your prompt answer.
> A related issue:
> I am interesting in getting a vector of points values from the SPDE

> Mat�rn model on S^2

> where the sampling sites belong to a small area on the sphere (large
> enough for the curvature to matter).
> Borrowing from http://www.r-inla.org/examples/case-studies/simpson2011
> I got something along the line of the below which seems to work fine
> for a window large enough
> but returns many NAs for a small window. Any work around?

As far as I can tell, those NAs are entirely correct and shouldn't be
worked around. Explanation below.
From looking at this code:

> ## proj =
> inla.mesh.projector(mesh,dims=c(180,90),projection="longsinlat")
> proj = inla.mesh.projector(mesh,dims=c(180,90),projection="longlat")
> ## better than longsinlat in terms of not getting NAs

> <cut>

> plotdata.x = inla.mesh.project(proj, x)
> mean(is.na(plotdata.x))

It appears that you're trying to project a function on the mesh onto a
planar projection (using longsinlat, which is an equal-area projection,
or longlat, which is not) for a lattice of points in the projection
domain. In all cases where the boundary of the mesh is not projected
exactly onto the borders of that lattice (i.e. in nearly _all_ cases,
since a typical mesh is a non-rectangular convex set) you will then get
NA for all points in the lattice which corresponds to a point outside
the mesh. Thus, the NA:s you get should be there, as the model is
completely undefined outside of the mesh. This is the intended
behaviour of inla.mesh.project(), as it interacts well with plotting
methods such as levelplot(), which will (correctly) ignore the NA "values".

if the above does not answer your question, I would need a more specific
explanation of what you are trying to accomplish, or if you are getting
NAs in places one shouldn't expect them, given the above.

Finn

[Finn Lindgren]

Hi again,

after running your code I think I see what the issue is; the function
that determines the automatic limits of the projection lattice thinks
that it needs the entire globe even though the mesh is really small:

inla.mesh.map.lim(loc = mesh$loc, projection = "longlat")
$xlim
[1] -180 180
$ylim
[1] -90 90

This is clearly not what you want in your case. I'll try to figure out
why those automatic values are so strange. In the meantime, if you
instead explicitly specify xlim and ylim in your call to
inla.mesh.projector() you will have direct control over which region the
lattice should cover.

Finn

[Finn Lindgren]

On 19/03/13 15:07, Finn Lindgren wrote:
>
> inla.mesh.map.lim(loc = mesh$loc, projection = "longlat")
> $xlim
> [1] -180 180
> $ylim
> [1] -90 90

...and now, having looked at what inla.mesh.map.lim() does, the
explanation is simple:
It makes no attempt whatsoever to find a small covering rectangle for
spherical meshes, but _always_ gives the whole globe.
Therefore you'll need to specify the projection limits explicitly.

The reason it does not attempt to do anything more clever is to avoid
surprises for meshes crossing the -180/180 degree longitude line. For
nice and small convex domains it would be possible, but defining a
sensible method that covers all non-convex and non-connected domains is
hard. The users themselves however will have a much better idea of what
their domain means and where it is located.

Finn