r - Create grid (or polygons) over shapefile that vary in area but equal in number of points encompassed

Question

I am using R's various spatial analysis features. I have data from a number of sampling locations throughout Africa, distributed unevenly. For several reasons, it is of interest to break our Africa map down into grids and compare characteristics of grids, matched on several characteristics. Ultimately, we want to be able to randomly sample grids. However, due to the uneven distribution of sampling locations, having evenly sized grid cells is problematic, because we would end up with some grids that have only one data point and other grids that have hundreds.

Basically, I know how to generate grid cells that are equal in area, with varying numbers of sampling locations ... however, what I want is essentially the opposite. Grid cells that vary in AREA, but with equal (or, realistically, near-equal, within a certain range) numbers of sampling locations. I recognize, too, that the solution may very likely involve polygons, rather than rectangular grid cells, but I am not sure how to implement this.

Simplified, reproducible example below, which takes a shapefile of Africa, randomly generates locations across Africa (ignoring for now that many of these are randomly generated in the ocean), and layers an even grid across it:

library(maptools)
library(lattice)

# Import shapefile of Africa
data(wrld_simpl, package="maptools")
africa <- wrld_simpl[wrld_simpl$REGION==2,]

# Randomly generate 200 coordinates
locations <- cbind.data.frame( x = runif(200, min(coordinates(africa)[,1]), max(coordinates(africa)[,1])),
                               y = runif(200, min(coordinates(africa)[,2]), max(coordinates(africa)[,2])))
coordinates(locations) <- ~ x + y
proj4string(locations) <- proj4string(africa)

# Verify with plot
plot(africa)
points(locations, pch=16, bg="black")

# Create spatial grid
bb <- bbox(africa)
cs <- c(10, 10)  # cell size 
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
sp_grd <- SpatialGridDataFrame(grd,
                               data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(africa)))

# Plot grid with randomly generated coordinates
spplot(sp_grd, "id", colorkey=FALSE,
       panel = function(..., col.regions) {
         panel.gridplot(..., border="black",col.regions="white")
         sp.polygons(africa)
         sp.points(locations, cex=1.2, pch=16, col="black")
         panel.text(...,col="red")
       })

# Can use 'over' funciton to find grid cells associated with each of our random points
# over(locations,sp_grd)

So, I know how to overlay this grid, count the number of sampling locations in each cell, extract data values for each cell, etc. etc. However, I am not sure how to proceed from here.

How would I generate grid cells (or polygons) that encompass the extent of my map area of interest, whereby the geographical area of the cells (or polygons) is variable but the number of sampling locations within each is fixed (or within a certain range)?

I am very new to the spatial analysis world.

Answer 1

Here's a method that first chops the data by Y coordinate forming horizontal bands of roughly equal numbers of points, and then chops each band into slices containing roughly equal numbers of points.

It takes a data fame with $x and $y columns, and N, the rough number of points you want in each cell:

chop <- function(xy,N){
    Npts = nrow(xy)
    Nyh = round(sqrt(Npts/N))
    Nxh = Nyh
    yp = seq(0,1,len=Nyh+1)
    ybreaks = quantile(xy$y,yp)
    ybins = cut(xy$y, ybreaks, include.lowest=TRUE, labels=FALSE)
    xy$ybin = ybins
    xy$ymin = ybreaks[ybins]
    xy$ymax = ybreaks[ybins+1]
    for(ybin in ybins){
        inbin = ybins==ybin
        xyb = xy[inbin,"x"]
        xp = seq(0,1,len=Nxh+1)
        xbreaks = quantile(xyb,xp)
        xbins = cut(xyb, xbreaks, include.lowest=TRUE, labels=FALSE)
        xy$xmin[inbin]=xbreaks[xbins]
        xy$xmax[inbin]=xbreaks[xbins+1]
        xy$xbin[inbin] = xbins
        }
    xy$bin = paste0(xy$xbin,",",xy$ybin)
    xy
}

# helper to get cell coordinates
cells <- function(chopped){
    xy = chopped[!duplicated(chopped$bin),c("xmin","ymin","xmax","ymax")]
    xy
}

And here's a test:

> set.seed(123)
> xy = data.frame(x=runif(123, 100,200), y=runif(123,200,400))
> plot(xy$x, xy$y, asp=1)
> grid = chop(xy, 10)
> head(grid)
         x        y ybin     ymin     ymax     xmin     xmax xbin bin
1 128.7578 243.9535    1 202.0934 253.0719 128.0981 165.3102    2 2,1
2 178.8305 273.8978    2 253.0719 296.4085 170.0444 195.4474    4 4,2
3 140.8977 396.8438    4 338.8269 397.1282 137.1654 154.4066    2 2,4
4 188.3017 230.8405    1 202.0934 253.0719 165.3102 189.5919    3 3,1
5 194.0467 218.2088    1 202.0934 253.0719 189.5919 199.4270    4 4,1
6 104.5556 228.3814    1 202.0934 253.0719 102.4614 128.0981    1 1,1

Here bin is the index of the resulting bin, with its coordinates in the xmin etc columns. xbin and ybin are just the parts of that.

You can see how many points are in each bin:

> table(grid$bin)

1,1 1,2 1,3 1,4 2,1 2,2 2,3 2,4 3,1 3,2 3,3 3,4 4,1 4,2 4,3 4,4 
  8   8   8   8   8   8   7   8   7   7   7   7   8   8   8   8 

The algorithm for bin numbers isn't perfect, so in this case you don't get 10 points per bin. This is because I'm using equal numbers of X and Y bins, so 123 points in 16 (irregular) bins is 7.68 per bin. A better method for choosing Nxh and Nyh could improve this (in this case, Nxh=4 and Nyh=3 would give you 10.25 points per bin).

The cells function takes the output and returns the coordinates for each of the grid cells (16 here) and you can plot them with rect:

> p = cells(grid)
> rect(p$xmin,p$ymin,p$xmax,p$ymax)
> 

Giving: