I have a SpatRaster where cells are either NA or 1. There are some points on it.

r <- rast(ncols=36, nrows=18, crs="+proj=longlat +datum=WGS84")
r[1:200] <- 1
plot(r)

p2 <- vect(rbind(c(30,-30), c(25,40), c(-9,-3)), crs="+proj=longlat +datum=WGS84")
plot(p2, add = T)

Is there a way to calculate the distance of only the non-NA cells to the points? I'm trying to reduce computation when calculating distances on a very large raster that is mostly NAs.

I tried the following:

dist<- distance(x = not.na(r), y = p2)
plot(dist)

This is my desired result:

I have a SpatRaster where cells are either NA or 1. There are some points on it.

r <- rast(ncols=36, nrows=18, crs="+proj=longlat +datum=WGS84")
r[1:200] <- 1
plot(r)

p2 <- vect(rbind(c(30,-30), c(25,40), c(-9,-3)), crs="+proj=longlat +datum=WGS84")
plot(p2, add = T)

Is there a way to calculate the distance of only the non-NA cells to the points? I'm trying to reduce computation when calculating distances on a very large raster that is mostly NAs.

I tried the following:

dist<- distance(x = not.na(r), y = p2)
plot(dist)

This is my desired result:

• r
• distance
• raster
• terra

• I am not quite sure what your gripe is with your own solution. Can you provide an example of the desired result? – D.J Commented 21 hours ago
• 2 Thanks, I've updated the question above. – patch105 Commented 21 hours ago

3 Answers 3

You could use the exclude parameter of distance if you use a single raster. To do that, set the cells that overlap the target points to a value different from the rest of the raster:

library(terra)
# create dataset
r <- rast(ncols=36, nrows=18, crs="+proj=longlat +datum=WGS84")
r[1:200] <- 1
p2 <- vect(rbind(c(30,-30), c(25,40), c(-9,-3)), crs="+proj=longlat +datum=WGS84")

# set values of cells overlapping points to -1
r[cells(r,p2)[,2]] <- -1
dist <- distance(x = r, target = 1, exclude= NA)
plot(dist)

Here is a solution using sf/stars as I am proficient with those but have little experience with the terra package.

library(terra)
library(sf)
library(stars)
library(ggplot2)

# Your Code
r <- rast(ncols=36, nrows=18, crs="+proj=longlat +datum=WGS84")  
r[1:200] <- 1  
p2 <- vect(rbind(c(30,-30), c(25,40), c(-9,-3)), crs="+proj=longlat +datum=WGS84")

# Convert to sf/stars
raster_stars <- st_as_stars(r)
points_sf <- st_as_sf(p2)  

raster_coords <- st_as_sf(as.data.frame(st_coordinates(raster_stars)), 
                          coords = c("x", "y"), 
                          crs = st_crs(raster_stars))

valid_coords <- raster_coords[which(!is.na(raster_stars[[1]])), ]  # Remove NA pixels

# Compute distance matrix
dist_matrix <- st_distance(points_sf, valid_coords$geometry)

min_distances <- apply(dist_matrix, 2, min)  # Column-wise min (distance to closest point)
valid_coords$distance <- min_distances  # Attach distances to raster coordinates

# Convert distance to raster 
distance_raster <- st_as_stars(st_rasterize(valid_coords, template = raster_stars))

# ggplot2 version with color scale and point overlay
x11(); ggplot() +
  geom_stars(data = distance_raster) +
  scale_fill_viridis_c(name = "Distance (m)") +  # Color scale
  geom_sf(data = points_sf, color = "red", size = 3) +  # Overlay points
  ggtitle("Distance from Points to Raster Pixels") +
  theme_minimal()

Calculating distance() using p2 and applying mask() afterwards seems to be faster than trying to exclude cells within the SpatRaster. It does not answer your "how to exclude cells question", but it does give useful context.

library(terra)

# Larger example of your SpatRaster
set.seed(42)
r <- rast(ncols = 3600, nrows = 1800, crs = "+proj=longlat +datum=WGS84")
r[1:2e6] <- 1

# Example points
p2 <- vect(rbind(c(30, -30), c(25, 40), c(-9, -3)), crs = "+proj=longlat +datum=WGS84")

# Caculate distance and mask()
dist <- distance(r, p2)
dist <- mask(dist, r)

plot(dist)

My (corrected) original answer involves converting r to a point SpatVector, however this was not faster than the above method:

# Convert r to points and calculate distance
v <- as.points(r)
dist <- distance(v, p2)

# Return minimum distance to any of the points
v$distance <- apply(dist, 1, min)

# Convert back to raster
dist <- rasterize(v, r, field = "distance")

Comparison:

library(microbenchmark)

microbenchmark(
  
  andrea = {
    r[cells(r, p2)[, 2]] <- -1
    dist <- distance(r, target = 1, exclude = NA)
  },
  
  lt1 = {
    dist <- distance(r, p2)
    dist <- mask(dist, r)
    
  },
  
  lt2 = {

    v <- as.points(r)
    dist <- distance(v, p2)
    v$distance <- apply(dist, 1, min)
    dist <- rasterize(v, r, field = "distance")
    
  },
  
  times = 10
)

# Unit: seconds                             
#    expr       min        lq      mean    median        uq       max neval cld
#  andrea 19.915640 20.102918 21.464082 20.598848 23.083761 25.725486    10 a  
#     lt1  2.139878  2.142092  2.212808  2.169019  2.219029  2.468196    10  b 
#     lt2 10.648845 11.274298 12.377071 12.239616 13.002600 14.802744    10   c