Count points within a polygon

Question

I'm trying to count the number of points within a polygon of a fixed size in R. I'm planning to do it a few times so that I can eventually calculate average density of points with some measure of error. The way I'm going about it is:

1. Create a 1 ha box
2. Have the box appear randomly within my study site
3. Count and save the number of points within the box
4. Rinse and repeat a few (hundred? thousand?) times
5. Calculate average density

My Google search turns up functions that generate random points within a polygon (e.g. st_sample, spsample), which isn't what I want. The closest I think I got was quadratresample but I couldn't figure out how to convert my SpatialPointsDataFrame to a ppp object. Even if I managed to do that, I'm not sure if I'm barking up the right tree with that function.

Is there a function that can count points from a polygon of a fixed size?

EDIT: I've posted and deleted another question about whether this process was the best way forward (after being advised that the post should be split into two separate questions). Jeffrey Evans has since killed two birds with one stone by offering a better route by way of PPA, as well as answering the specific question I had.

Answer 1

You are referring to a simple type of point pattern analysis (PPA) called Quadrat analysis. Quadrat analysis can be used both for assessing goodness of fit (of empirical data to a theoretical model) and for assessing the independence of two distributions. A chi-squared statistic is commonly used in evaluating the hypothesis of a random point process, using a Poisson distribution as the spatial null. However, there are other test statistics such as Person's X2 and Cressie-Read.

This analysis can be done in the spatstat package but, functions such as spatstat::quadratcount uses a normal gridded set of quadrats. With a custom function, it is "possible" to create random quadrats but, it really is not necessary because the quadrat.test function can perform a simulation, thus displacing the need for randomized quadrats because, randomization is occurring within the point process via a Monte Carlo simulation implemented. This is much more robust than what you are thinking. Traditional random quadrat tests are for testing specific types of hypotheses of spatial process. However, you may want to test a range of different sized quadrats to see what effect the quadrat sampling scheme is having on the results. This can have a notable effect on the results and is a know limitation of this approach.

If you use the maptools package coercing an sp object to ppp is as easy as as(x, "ppp") however, I will illustrate the long way around. Here is an example quadrat analysis using the meuse data set.

Add packages and data

library(sp)
library(spatstat)

data(meuse)
  coordinates(meuse) <- ~x+y

Create a convex hull window for the point pattern and coerce to a spatstat ppp object.

win <- spatstat::convexhull.xy(sp::coordinates(meuse))
meuse <- ppp(coordinates(meuse)[,1], coordinates(meuse)[,2], win)

Here we look at some basic statistics such as count and intensity in a given quadrat sampling scheme.

# counts 
( Q <- quadratcount(meuse, nx= 6, ny=6) )
 plot(meuse, pch=20, cols="grey70", main="Quadrat counts")
   plot(Q, add=TRUE)

# intensity
( Qd <- intensity(Q) )
plot(intensity(Q, image=TRUE), main=NULL, las=1)  
  plot(meuse, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)  

Now, we can test for dispersion using Complete Spatial Randomness (CSR) as our null. We will use a Monte Carlo approach with Cressie-Read power divergence as our test statistic.

( qtest <- quadrat.test(meuse, nx= 6, ny=6, method="MonteCarlo", nsim=999) )

Now, to directly address your question of randomly creating a box and counting points within. You can easily wrap this all up in a for loop to "rinse and repeat".

library(sp)
library(raster)
library(spatstat)
library(rgeos)
data(meuse)
  coordinates(meuse) <- ~x+y

For your sampling space, I am assuming that you want a convex hull. Now, here is where we create a random quadrats. The trick I am using here is to buffer a random point sample to the distance that would equal the desired size of the quadrat and then turn the extent of the buffer into a polygon. This gives me a square buffer. We do this n times, keeping track of counts, and store in a list object. The polygons can then be combined using do.call. Personally, I would randomly sample a range of quadrat sizes.

e <- gConvexHull(meuse)
s=250  # radius defining fixed size quadrat
# s = seq(100,300,10) # radius defining range sized quadrats
n=500  # 3 number of permutations
quadrats <- list()
  for(i in 1:n) {
    ss = sample(s,1)
    p <- as(extent(gBuffer(spsample(e, 1, "random"),width=ss)), "SpatialPolygons")
      quadrats[[i]] <- SpatialPolygonsDataFrame(p, data.frame(ID=i, 
                         n = length(which(meuse %over% p == TRUE))))
  }
quadrats <- do.call("rbind", quadrats)
  spplot(quadrats, "n")

If you want to add some rotation into the mix you can use the elide function in maptools.

p <- as(extent(gBuffer(spsample(e, 1, "random"),width=250)), "SpatialPolygons")
  p45 <- maptools::elide(p, rotate=45, center=apply(bbox(p), 1, mean))
    plot(p45)
      plot(p,add=TRUE)