I have a 2-dimmensional unit grid, and a bunch of line segments that start and end at any rational number. I need an efficient way to calculate which grid cells the line passes through. For example, the line:
From (2.1, 3.9) to (3.8, 4.8) passes through the grid cells with lower left points (2, 3), (2, 4), and (3, 4).
Is there a quick, efficient way to calculate these quadrants from the line's endpoints?
I'll be working in R, but an answer in Python or pseudocode would work too. Thanks!
Folks who work with spatial data deal with this kind of question all the time, so it may be worth piggy-backing on their efforts. Here's a solution that uses R's raster package (and functions from the sp package on which it depends):
library(raster)
## Create a SpatialLines object
a <- c(2.1, 3.9)
b <- c(3.8, 4.8)
## Method #1 -- Uses functions from the sp package.
SL <- SpatialLines(list(Lines(list(Line(rbind(a,b))), "ab")))
## Method #2 -- Uses readWKT() from the rgeos package. Easier to read.
# library(rgeos)
# string <- paste0("LINESTRING(", paste(a, b, collapse=", "), ")")
# SL <- readWKT(string)
## Create a raster object
m <- 10
n <- 10
mat <- matrix(seq_len(m*n), nrow = m, ncol = n)
r <- raster(mat, xmn = 0, xmx = n, ymn = 0, ymx = m)
## Find which cells are intersected & get coordinates of their lower-left corners
ii <- extract(r, SL, cellnumbers=TRUE)[[1]][, "cell"]
floor(xyFromCell(r, ii))
# x y
# [1,] 2 4
# [2,] 3 4
# [3,] 2 3
## Confirm that this is correct with a plot
image(r)
plot(as(rasterize(SL, r), "SpatialPolygons"),
border = "darkgrey", lwd = 2, add = TRUE)
lines(SL)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With