Check if a point lies on the boundary of two raster cells-CodePudding

I have 25-km raster grid and a point. If I plot the point on the raster, it falls ~on the borderof two adjacent raster cells. I don't how to generate a sample data for this specific case but here's a visualisation of what I am taking about.

Is there any way I can detect for a group of points how many points are cases like this which could like very close to the border of two adjacent raster cells?

CodePudding user response：

You have the raster extent, dimensions and resolution, you can calculate the coordinates of the cells, like

library(terra)
#> terra 1.6.17
r <- rast(nrows=20, ncols=10, xmin=0, xmax=10, ymin = 0, ymax = 10)
for (i in seq_len(nrow(r))) {
  for (j in seq_len(ncol(r))) {
    if(i %% 2 != 0 & j %% 2 != 0) {
      print(paste(i, j))
      r[i,j] <- 1
    }
  }
}
for (i in seq_len(nrow(r))) {
  for (j in seq_len(ncol(r))) {
    if(i %% 2 == 0 & j %% 2 == 0) {
      print(paste(i, j))
      r[i,j] <- 1
    }
  }
}

terra::plot(r, axes = TRUE)

# a few points
p <- rbind(c(1,1), c(2.3,3), c(3, 4), c(5.5,5.5)) |>
  terra::vect("points")
terra::plot(p, add = TRUE, pch = 12)

Now the calculations. Dimension and resolution:

dim(r)
#> [1] 20 10  1
res(r)
#> [1] 1.0 0.5

Let's calculate the coordinates of the cells (here just X):

xcoords <- vector(mode="list", length=dim(r)[1]-1)
for (i in seq_along(xcoords)) {
 xcoords[[i]] <- as.numeric(ext(r)[1]) i*res(r)[1]
}

And finally chceck, which points are located exactly on the xcoords:

xx <- geom(p)[,"x"]

which(xx %in% xcoords)
#> [1] 1 3

Having raster dimensions you can create lines corresponding to grid, and do spatial intersects, see sf::st_make_grid().

^{Created on 2022-11-20 with}

To get the distance to the nearest border

to_border <- function(r, x, y) {
    i <- (x >= xmin(r)) & (x <= xmax(r)) & (y >= ymin(r)) & (y <= ymax(r))
    d <- rep(NA, length(i))
    d[i] <- (x[i] - xmin(r)) %% res(r)[1]
    d[i] <- pmin(d[i], (y[i] - ymin(r)) %% res(r)[2])
    d
}

to_border(r, x, y)
#[1] 0.0 0.5 0.0 0.0  NA

In principle, you could also create lines from the raster, and do a relate query like I show below, but that is really inefficient.

pts <- vect(cbind(x, y), crs=crs(r))
rlns <- aggregate(as.lines(r))
relate(pts, rlns, "intersects")
#      [,1]
#[1,]  TRUE
#[2,] FALSE
#[3,]  TRUE
#[4,]  TRUE
#[5,] FALSE

But if you wanted the distance to the border, you may need to do something like this.

distance(pts, rlns)
#     [,1]
#[1,]  0.0
#[2,]  0.5
#[3,]  0.0
#[4,]  0.0
#[5,]  0.2