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

Microsporophyll

Maybe first (counter)-rotate the polygon, then cover it with a grid, which you then rotate to cover the original polygon? – Uvarovite

good idea! I'll give it a try. – Microsporophyll

Unfortunately, @JoshO'Brien suggestion did not work for me... still getting incomplete cover... Anyone else? – Microsporophyll

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)

Henotheism

Indeed you are right! I was rotating "back" the grid using the centroid of the rotated polygon instead than the centroid of the "original" one. Your solution seems to work perfectly. Thank you very much! – Microsporophyll