Creating a regular polygon grid over a spatial extent, rotated by a given angle

Question

Hi all,

I am struggling with this and hope someone could come out with a simple solution.

My objective is to create a regular polygon grid over the extent of a polygon, but rotated by a user-defined angle.

I know that I can easily create a North/South polygon grid in sf using for example:

library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3
inpoly <- st_read(system.file("shape/nc.shp", package="sf"))[1,] %>% 
  sf::st_transform(3857) %>% 
  sf::st_geometry()
grd <- sf::st_make_grid(inpoly, cellsize = 3000)
plot(inpoly, col = "blue")
plot(grd, add = TRUE)

I also know that I can easily rotate it by a given angle using:

rotang = 20
rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
grd_rot <- (grd - st_centroid(st_union(grd))) * rot(rotang * pi / 180) +
  st_centroid(st_union(grd))
plot(inpoly, col = "blue")
plot(grd_rot, add = TRUE)

My problem is that , depending on the rotation angle, the general “orientation” of the input polygon and the cell size, the rotated grid may not cover anymore the full extent of the polygon, as shown below:

rotang = 45
rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
grd_rot <- (grd - st_centroid(st_union(grd))) * rot(rotang * pi / 180) +
  st_centroid(st_union(grd))
plot(inpoly, col = "blue")
plot(grd_rot, add = TRUE)

Any clever idea about how I could address this issue and create a rotated grid fully covering the polygon (besides by creating a larger grid to start with, which is quite inefficient for small cellsizes?)?

Either sf or sp solutions would be welcome. “Bonus points” if it is possible to make the grid start at one of the extreme vertexes of the polygon (i.e., the first line of the grid “touches” the northern vertex of the polygon), but that is not "mandatory".

Created on 2018-07-11 by the reprex package (v0.2.0).

Answer 1

You didn't specify, how exactly @JoshO'Brien's suggestion doesn't work for you, but I suspect you rotated polygon and grid around different rotation centers. You didn't specify any constraints on rotation origin, so I assume in my code snippet below that it's not important, but you can use any point as soon as it's the same for both rotations:

library(sf)
rotang = 45
rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
tran = function(geo, ang, center) (geo - center) * rot(ang * pi / 180) + center
inpoly <- st_read(system.file("shape/nc.shp", package="sf"))[1,] %>% 
  sf::st_transform(3857) %>% 
  sf::st_geometry()
center <- st_centroid(st_union(inpoly))
grd <- sf::st_make_grid(tran(inpoly, -rotang, center), cellsize = 3000)
grd_rot <- tran(grd, rotang, center)
plot(inpoly, col = "blue")
plot(grd_rot, add = TRUE)