Creating hexagons around polygon using sf

Question

I create hexagons within a polygon perimeter using either the sp package or the newer sf package. Page 134 6.4 of Bivand et al.'s (2008) "Applied Spatial Data Analysis with R" gives an example of creating hexagons with sp that is repeated in form here:

StudyArea <- read_sf(paste0(here("Inputs"),"/StudyArea.shp"))
StudyArea_sp <- sf::as_Spatial(StudyArea $geometry)

HexPts <- spsample(StudyArea_sp, type = "hexagonal", cellsize = 21750)
HexPols <- HexPoints2SpatialPolygons(HexPts)
HexPols <- sf::st_as_sf(HexPols)

This creates the following output:

I get a similar, but not the exact same result with sf style:

HexPols2 <- st_make_grid(StudyArea , cellsize = 21750, square = FALSE)
HexPols3 <- st_intersection(HexPols2,StudyArea)
plot(StudyArea)
plot(HexPols3, add=TRUE)

I have two interrelated questions:

1. There is a difference between the sp and sf approach that results in a different fitting of the hexagons. Why would this occur and how should I best control this? For example, I would like to maximize the number of full hexagons that occur within the polygon and to minimize the number of hexagons that touch the perimeter.

2. How do I get the hexagons to give complete coverage of the polygon? The sp approach has complete hexagons, but it does not always touch the perimeter. In contrast, the sf approach touches the perimeter but does not finish drawing the hexagons that meet the perimeter. I would like hexagons that touch the perimeter to be completed and not cut off at the boundary. However, I have tried several options, such as st_overlaps, and get the following type of error:

   Error in xy.coords(x, y, xlabel, ylabel, log) : 'x' is a list, but does not have components 'x' and 'y'

Answer 1

For Q1 (a) (why is sp/sf different?) the difference is just down to however the algorithm decides to start its hexagonal grid (either in points for spsample or polys for st_make_grid). For Q1 (b) (maximising hexes within polygon) you'd create a large hexagonal grid that spans right over the polygon and beyond, then count the number of hexagons within the polygon for shifts in X and Y over the lattice size of the hexagons. This is then a 2d maximisation problem, but its non linear and there's no easy way to find the global maximum. Also note we like to stick to one question per post, so this should maybe be cut out of here and made into a new question if you want a full answer.

For Q2, in sp: you could buffer the polygon, generate the points with spsample and then run HexPoints2SpatialPolygons (where does that come from?) on this extended hexagonal lattice point set, then crop to the polygons. You'd need the rgeos package for all this, but its probably not worth the bother with the sp package etc being retired shortly. You should maybe cut this out of the question too.

For Q2 in sf: generate the polygons as you do, but use st_intersects to test if the hexagon is over any part of the polygon. Then subset.

# get sample data from package
nc = st_read(system.file("shape/nc.shp", package="sf"))
# convert from lat-long to a cartesian - use a better CRS for your data
nc = st_transform(nc, 3857)

# make a hex grid
g = st_make_grid(nc, square=FALSE)

# test if hex grid intersects any polygons. st_intersects
# returns a list, and the element list is zero if there's no intersection
ig = lengths(st_intersects(g, nc)) > 0

# plot the map, and add the intersected hexagons:
plot(st_geometry(nc))
plot(g[ig], col="red", add=TRUE)