Creating a rotated but regular points grid (with r)

Recently at work we have had to deal with a situation that has given us a little headache.

A normal common situation (at least for us) is to have a matrix or spatial layer of points in a “regular” situation (let’s call it that way), or in a “straight” form (orthogonal) to the reference axes.

However, for the sake of life (read, regular mesh for modeling coastal dynamics), we needed a regular rotated mesh. The important term (the really important one) is rotated. The “direct” option (the first thing that comes to mind) is to create a straight regular mesh, and turn it. Solved. but the reality is that it is not something so simple. Rotation is one type of “affine transformations”. An affine transformation, which includes among others, shifting (translation), scaling and rotation, is any transformation that preserves lines and parallelism. However, angles or length are not necessarily preserved (Lovelace, 2019). And trust me, lengths are usually not preserved.

The situation (and I don’t want to go into more detail about the contour conditions) was also a bit more complicated because we necessarily had to work in UTM coordinates. So in short, what we wanted and needed was a regular mesh, but rotated with respect to the coordinate axes.

Let’s get to work. If we can’t create a regular mesh and then rotate it, let’s use the power of trigonometry. Our question was: are we capable of knowing the coordinates of a point from a given origin, knowing the dimensions of the mesh, and the angle of rotation? The trigonometric answer: yes. So after some comings and goings, we arrive at the code that is attached below.

upper_left_x <- 5000
upper_left_y <- 5000
L <- 2000 # dimension from upper_left to bottom_left
dx <- 200 # delta in L dimension
M <- 1000 # dimension from upper_left to upper_right
dy <- 200 # delta in M dimension
angle <- pi/6  # angle between M dimension and horizontal. IMPORTANT IN RADIANS


coords_points <- data.frame(x_point = double(),
                            y_point = double(),
                            stringsAsFactors=FALSE)

# Delta L
deltaL <- coords_points
for (i in seq(0, L, dx)){
  deltaL_x_point <- + dx*i*cos(pi/2 -angle)
  deltaL_y_point <- - dx*i*sin(pi/2 - angle)
  
  deltaL <- rbind(deltaL, c(deltaL_x_point, deltaL_y_point))
}

# Delta M
deltaM <- coords_points
for (i in seq(0, M, dy)){
  deltaM_x_point <- dy*i*cos(angle)
  deltaM_y_point <- dy*i*sin(angle)
  
  deltaM <- rbind(deltaM, c(deltaM_x_point, deltaM_y_point))
}


coords_dataframe <- as.data.frame(cbind(rep(upper_left_x, (L/dx+1)*(M/dy+1)), rep(upper_left_y, M/dy+1)))
deltaL_dataframe <- deltaL[rep(seq_len(nrow(deltaL)), times=M/dy+1),]
deltaM_dataframe <- deltaM[rep(seq_len(nrow(deltaM)), each=L/dx+1),]


coords_points <- coords_dataframe + deltaL_dataframe + deltaM_dataframe


colnames(coords_points) <- c("x_point", "y_point")
# OUTPUT 
write.csv(coords_points, file = "output_path/file.csv") # Define output path

The most important thing about it is: IT WORKS.

library(sp)
spdf <- SpatialPointsDataFrame(coords = coords_points, data = coords_points,
                               proj4string = CRS("+init=EPSG:32618")) # Define own reference system
plot(spdf)

Obviously we started by creating a code based on for cycles. It worked too, but its execution time was quite slow (more than 24 hours later it had only executed 30%). So we looked for an alternative way. The execution time has now dropped to just minutes, with the export process taking the longest part of those minutes.

About the code itself:

• First lines defines the geometry of the problem.
• An auxiliary data frame is then created.
• Deltas matrices are calculated. These includes the differences in x and y coordinates from the origin point.
• The final coordinates of all points are obtained, starting from origin point and summing deltas matrices.
• Wrapping up, saving and verifying output.

As we are basically good people, we wanted to share it here, in case someone in some moment of desperation finds it and can save a life.