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

Maybe first (counter)-rotate the polygon, then cover **it** with a grid, which you then rotate to cover the original polygon? — Josh O'Brien, Jul 11 '18 at 13:37
Unfortunately, @JoshO'Brien suggestion did not work for me... still getting incomplete cover... Anyone else? — lbusett, Jul 16 '18 at 16:42

score 12 · Accepted Answer · answered Jul 17 '18 at 05:24

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)

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! — lbusett, Jul 17 '18 at 08:37

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

1 Answers1

Linked