sample.lgcp() returns matrix, not SpatialPoints object

[david-borchers]

library(inlabru) library(INLA) library(ggplot2)

(First Unzip Archive.zip, getting the files ingulaboundary.rds and meshspdf.rds)

Archive.zip ingulaboundary <- readRDS("ingulaboundary.rds") # boundary file meshspdf <- readRDS("meshspdf.rds") # SpatialPointsDataFrame giving locations to use for mesh plot(ingulaboundary); plot(meshspdf,add=TRUE) # plot to see what we have mesh0 <- inla.mesh.2d(loc=meshspdf, boundary=ingulaboundary, max.n=500) ggplot() + gg(mesh0) + gg(ingulaboundary) + coord_fixed(ratio = 1) # plot it to check

Draw a sample of points from a Poisson process with log intensity given by loglambda vector:

set.seed(1) # for reproducibility loglambda = log(rep(1e5, mesh0$n)) popts = sample.lgcp(mesh0, loglambda=loglambda) class(popts)

[finnlindgren]

The code currently only returns a SpatialPoints object if the CRS is set on the mesh. Since SpatialPoints CRS info can be NA and the documentation explicitly claims it should return SpatialPoints, I'll change it to do that. I'll check there isn't documentation examples relying on matrix output as well, since indexing doesn't work the same for matrices. A simple as.data.frame (but not as.matrix) can be used to work around such use.

[finnlindgren]

The code behaves differently for different strategy settings; for rectangle, any mesh coordinate system is ignored, and the "raw" 2D space is used. For spherical and spliced-spherical, it samples point on a globe, assuming that lambda is given in counts per square km and transforms them to longlat (and then onward to the mesh crs). If the mesh doesn't have a crs, and the mesh isn't a sphere, since then the first transformation must fail. Also, "sliced-spherical" assumes the input is in longlat when finding the coordinate range, but doesn't check that.

What I think it should do:

• For crs-less meshes on R2: use raw units and set NA crs on output
• For crs meshes on R2: interpret lambda as per km^2 (Allow the R parameter to override the default Earth radius) and set same crs on output
• For crs-less meshes on S2: interpret lambda with raw units. Allow the R parameter to override the sphere radius of the mesh. Output in Euclidean sphere coordinates, set geocentric crs with the same radius as the mesh (or R, if specified).
• For crs meshes on S2: interpret lambda in counts per km^2, Allow the R parameter to override the default Earth radius. Output in Euclidean sphere coordinates, and set geocentric crs with the same radius as the mesh (which should be the original mesh crs).
• Add a parameter ignore.CRS to override any mesh crs information (making the first and third version the only possible ones).

This means that the sampled points will live on the same mesh as the input mesh, whatever that is.

[david-borchers]

Thanks for the instant response Finn.

David