2

I have a SpatVect consisting of points and I want to rasterize them into a SpatRaster with a given resolution. Is there a way of specifying a function taking in the points that are within a buffer of each raster cell?

Many thanks Joao

-- Update -- Maybe a figure would help understand what I'm after with my question. The red square will have to be run over the center of each pixel to calculate some statistics using the ovelaying points. Apologies for the clumsy question, but I hope the figure is clear enough...

enter image description here

Joao Carreiras
  • 45
  • 5
  • The picture certainly helps. And what would the derived values be, ie. what statistics are sought? – Chris Sep 22 '22 at 21:45
  • First, just the number of points (length) but afterwards I will need to make more complex calculations, but for that I already have a function – Joao Carreiras Sep 22 '22 at 21:53
  • So multiple counts of the same point as the 3x3 window moves along doesn't matter. This is nearly `focal`, but `focal` doesn't use square windows, rather rectangular, but does allow a virtual raster approach (expand = TRUE) for when the rect box will be outside edges. I'll head scratch some, and in the meantime, someone far better informed will come by. Actually inclined to take this toward `sf` for some of the approaches `st_within` & etc, but head scratch there as well. – Chris Sep 22 '22 at 22:15

2 Answers2

2

terra version 1.6-28 supports rasterization of points with a rectangular moving window.

Example data

library(terra)
#terra 1.6.33
r <- rast(ncol=100, nrow=100, crs="local", xmin=0, xmax=50, ymin=0, ymax=50)
set.seed(100)
x <- runif(50, 5, 45)
y <- runif(50, 5, 45)
z <- sample(50)
v <- vect(data.frame(x,y,z), geom=c("x", "y"))

Solution

r1 <- rasterizeWin(v, r, field="z", fun="count", pars=10, win="rectangle")

plot(r1)
points(x, y)

enter image description here

You can change fun to another function that works for you, and you can change the size of the moving window with pars.

Instead of a rectangle, you can also use a circle or an ellipse. The border of a circular window is equidistant from the center of the cells. In contrast, the border of rectangles are at a constant distance from the border of the grid cells in most directions (not at the corners). Here is an example.

r2 <- rasterizeWin(v, r, field="z", fun="count", pars=5.25, win="circle")
plot(r2)

enter image description here

You can also use buffers around each cell to get a window that is truly equidistant from each cell border.

r3 <- rasterizeWin(v, r, field="z", fun=length, pars=5, win="buf")
plot(r3)

enter image description here

In this case, because the buffer size is large relative to the cell size, the result is very similar to what you get when using a circular window. Using "circle" should be the fastest, and using "buffer" should be the slowest in most cases. The function should now in all cases be memory-safe, except, perhaps when using very large buffers (more could be done if need be).

Version 1.6-28 is currently the development version. You can install it with install.packages('terra', repos='https://rspatial.r-universe.dev')

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Thank you Robert, this is perfect and just what I needed. – Joao Carreiras Oct 04 '22 at 08:48
  • I'm getting: Error in x@ptr$winrect(y[, 1], y[, 2], pars, opt) : std::bad_alloc Memory issues? Raster is ~ 8000x8000 pixels and vector ~3.5 M onservations. I'm runnign this on a machine with 195 GB RAM, of which 145 GB is still free. – Joao Carreiras Oct 04 '22 at 09:58
  • I have improved the memory safety and updated the answer. Can you please try again? – Robert Hijmans Oct 04 '22 at 20:34
  • Hi Robert. I'm using terra 1.6-28 and I'm still getting the same error message: Error in rbe@ptr$winrect(y[, 1], y[, 2], pars, opt) : std::bad_alloc Quick question: pars should be given in projected units or degrees? – Joao Carreiras Oct 05 '22 at 10:03
  • In the units of the CRS, so if your data is lonlat so should the pars be. Can you contact me privately and share the data and script so that I can have a look? – Robert Hijmans Oct 05 '22 at 15:37
  • Hi Robert. What happens at the border of the raster? Only the part of the window covering the raster will be used? – Joao Carreiras Oct 16 '22 at 14:53
  • The entire window is used. – Robert Hijmans Oct 20 '22 at 03:53
0

The approach you take seems to depend on what result you're looking for from the above and the relationship they have with each other.

library(terra)

`terra::buffer(` # both SpatVectx/SpatRastery, to distance in 'm'
`terra::buffer(` # that is meaningful
#take Rasty to SpatVecty 
`terra::as.polygons(`, #then 
`z<-terra::intersection(SpatVectx, SpatVecty)`

then back to SpatRastz? terra::mask or crop, might also be useful, again depending on where things are going next.

Chris
  • 1,647
  • 1
  • 18
  • 25