Documentation request: cellsize argument st_make_grid()

[JosiahParry]

Hi there!

I'm working to make a hexagonal tessellation using st_make_grid(). Perhaps this is an oversight on my end, but the argument cellsize is unclear to me. When specifying cellsize what am I specifying? Am I specifying the radius from the centroid? The width of a rectangle? Is the cellsize different depending on the square argument? e.g. if square = TRUE it uses rectangle width and if square = FALSE it uses circle radius?

If you can specify in this issue, I'd be happy to make a pull request to supplement the documentation.

Thanks!

[edzer]

This is a bit of guess work for me:

• cellsize[1] / sqrt(3) is passed as dx to make_hex_grid(); I don't know why, maybe because the area of the hex cell is similar or equal to that of a square cell with size cellsize[1]?
• that function mentions: "has x spacing dx: the shortest distance between x coordinates with identical y coordinate"
• I assume this is valid for flat_topped = TRUE, and refers to the case with x/y swapped when FALSE

Can you work it out from there?

[JosiahParry]

Thanks for pointing that out @edzer. That's actually rather helpful. This is a bit out of my expertise but I'll make my best guess work.

Following this SO thread's answer and what I see in the populartimes python library (which uses the honeycomb pattern for grid making) the sqrt(3) * r is the distance between the centroids of hexagons. The dx referenced in make_hex_grid() I suspect is supposed to be the width of the cell (square or hexagon). It is then divided by two to calculate the radius.

That's a long way of saying: my conclusion here is that it's the diameter of the square or rectangle.

Can someone else confirm?

[JosiahParry]

Proposed:

@param cellsize The diameter of the square or hexagon used to create a grid in the same units as the provided object.

[edzer]

We can also measure:

library(sf)
b = st_as_sfc(st_bbox(c(xmin=0,xmax=1,ymin=0,ymax=1)))
h = st_make_grid(b, cellsize = .1, square = FALSE)
plot(b)
plot(h, add = TRUE)

This shows that 0.1 is the shortest distance between two opposite ribs. Diameter is a property of a circle, for a square or hexagon I'm not sure what it refers to.

[kylebutts]

The "Diameter" of a regular polygon is technically two times the apothem: https://en.wikipedia.org/wiki/Apothem. But I don't think that's helpful in the context of the documentation. Perhaps "width" of the square or the hexagon will be clear?

[JosiahParry]

I think diameter would be the better parameter name. That could conceivable introduce breaking changes, however. I suspect my suggestion at https://github.com/r-spatial/sf/issues/1505#issuecomment-705718407 would suffice?

[edzer]

Yes, and in @details we should explain what we mean by diameter.

[rCarto]

In cartography::getGridLayer() I use a desired output area so the user enter something a bit more understandable for hexagons than diameter or width.

I don't know if it's relevant to add a new cellarea parameter to st_make_grid() though.

library(sf)
#> Linking to GEOS 3.7.1, GDAL 3.1.2, PROJ 7.1.0
library(units)
#> udunits system database from /usr/share/xml/udunits
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
nc <- st_transform(nc, "EPSG:6542")
# expected area = 10000m * 10000m = 1e+8 m2 = 100 km2
cellarea <- 100 * (1e+6)
cellsize <- 2 * sqrt(cellarea/((3*sqrt(3)/2))) * sqrt(3)/2
hexa <- st_make_grid(x = nc, cellsize = cellsize, square = FALSE)
a <- st_area(hexa)[1]
units(a) <- "km^2"
a
#> 100 [km^2]

^{Created on 2020-11-04 by the reprex package (v0.3.0)}

[edzer]

I like that! One would then require that either cellarea or cellsize is specified.

[plnnr]

This seems potentially useful:

[edzer]

I tried to improve the docs, and added support for area units to specify cellsize:

library(sf)
# Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
# hex:
b = st_as_sfc(st_bbox(c(xmin=0,xmax=1,ymin=0,ymax=1), crs = 3857))
h = st_make_grid(b, cellsize = .1, square = FALSE)
a = st_area(h)[1] # has units: m^2
h1 = st_make_grid(b, cellsize = a, square = FALSE)
all.equal(h, h1)
# [1] TRUE
# square:
nc <- read_sf(system.file("shape/nc.shp", package = "sf"))
nc <- st_transform(nc, 3857)
m1 = st_make_grid(nc, 100000)
m2 = st_make_grid(nc, units::as_units(100, "km")) # length units
all.equal(m1, m2)
# [1] TRUE
m3 = st_make_grid(nc, units::as_units(1e4, "km^2")) # area units
all.equal(m2, m3)
# [1] TRUE