Count points within a polygon

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)