# Creating hexagons around polygon using sf

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)
`````` 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'`

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