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 = ?))

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


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?

  • It is unclear to me what you want to achieve. What do you mean by "windows"? I do not understand what is the role of the polygon in this story. Can you draw a picture of what you want to to?
    – JRR
    Commented Jun 9, 2022 at 13:39
  • If I understand what you are wanting to do, calculate a PPA on a lidar point cloud given each pixel, this is an incorrect application of a point pattern statistics. In a planar context the statistic would be fairly irrelevant in indicating a spatial process and with the 3D nature of lidar data any PPA would be completely invalid (the CSR null does not hold). Commented Jun 9, 2022 at 15:46
  • @JRR In point-pattern analysis, the "window" is the area over which the process generating the point pattern is defined. ie its anywhere you might find a point from a realisation of the process. eg if you have a record of all trees in a field, the window is the shape of the field.
    – Spacedman
    Commented Jun 9, 2022 at 17:25
  • Why do you want to put the points into a regular grid? If you are interested in if your points are clustered then use a test from spatstat that works on the point coordinates. I don't understand "according to" in "I want to know the spatial distribution (sparse/clustered) according to a regular grid"
    – Spacedman
    Commented Jun 9, 2022 at 17:27
  • @Spacedman I know very little about point-pattern but here what I don't understand is the distinction between raster cells size and "windows". Does OP want to compute point pattern analysis of points inside each pixel with a windows smaller than the pixel size?
    – JRR
    Commented Jun 9, 2022 at 17:31

1 Answer 1


Use these packages:


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

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

enter image description here

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:

enter image description here

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:

enter image description here

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

enter image description here

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...)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.