Rasterise a point pattern in R

Question

I have a set of georeferenced points within an area for which I want to know the spatial distribution (sparse/clustered) according to a regular grid. I thought to describe the pattern using the Clark-Evans index (from spatstat.core package) so that I get a single value describing it...and this could be spatialised within the whole study area through a raster (I know that aggregation is scale-dependent but that's also part of the research interest).

Unfortunately, I can't achieve the rasterization process in a "smooth" way with R. The two approaches which I though of look something like:

my_fun = function(x,y){
         out = spatstat.core::clarkevans(spatstat.geom::ppp(x,y,window = ?))
         return(out)
         }


terra::rasterize(data, myraster, fun=my_fun)

or

lidR::pixel_metrics(data, fun=my_fun(X,Y), res=1)

Both approaches lack of a fondamental parameter: the window...and indeed points are dropped. ppp objects need a owin object which, in this case, it would be constituted by the pixel boundaries. But how can I provide a window at that stage? I guess creating a method?

Right now I'm using a pretty time-consuming workaround such as:

• convert raster pixels to polygons
• loop the aggregation function through each polygon/window
• revert back to raster making use of X,Y from the centroid and the index value I get.

Any smooth solution or methodological alternative?

Answer 1

Use these packages:

library(spatstat)
library(raster)
library(maptools)

create a test data set. 1000 points in a unit square and a 5x4 grid of quadrats:

nx=5
ny=4
grid = raster(matrix(1:(nx*ny), ny, nx))
gridpoly = rasterToPolygons(grid)
pts = cbind(runif(1000), runif(1000))
plot(gridpoly)
points(pts)

Now loop over each of the gridpoly rows, calling the spatstat function by creating a ppp with that row as the window, and it will warn (but drop) the points outside the window. Add the CE coefficient from the cdf corrected version to the polygons:

gridpoly$CE = sapply(1:nrow(gridpoly), function(i){
    clarkevans(as.ppp(pts, W=as.owin(gridpoly[i,])))["cdf"]
})

and plot:

giving a map of the coefficient for the points within each grid cell.

Testing with another data set that has a grid (extreme inhibitory) of points in most of the space and some random points in one corner:

gives a coefficient map which looks right - high values for inhibitory process except in the corner which is 1 for random points:

The difference in the blue area and slight difference in the yellow area is probably down to edge correction methods (see docs for the details).

It is be possible to do all this within spatstat using a split PPP object. For example here's a PPP of random points in a unit square:

P = ppp(pts[,1],pts[,2], window=owin(c(0,1), c(0,1)))

and this defines a grid of 5x4 in that window:

Z = quadrats(P, 5,4)

this Z can then be used to split the P object into parts which some spatstat function will apply over, for example density:

D = density(split(P,Z))

when plotted will show 20 separate density plots in a 5x4 layout. However a split ppp is also equivalent to a list, so you can apply any function over the splits:

CE2 = lapply(split(P,Z), clarkevans)

which gets you a list of clarkevans results:

$`Tile row 1, col 1`
    naive  Donnelly       cdf 
0.8908877 0.8346257 0.8620706 

$`Tile row 1, col 2`
   naive Donnelly      cdf 
1.277086 1.204744 1.217105 

etc which you can now manipulate with the usual list manipulation tools and reform into a grid. In fact, extracting the 3rd element from each of these gives the same vector as gridpoly$CE above:

 sapply(CE2, function(c){c[3]})
Tile row 1, col 1.cdf Tile row 1, col 2.cdf Tile row 1, col 3.cdf 
            0.8620706             1.2171046             1.1047117 

> gridpoly$CE
 [1] 0.8620706 1.2171046 1.1047117 0.864[etc]

which is reassuring that I've got both methods right (or both wrong...)