CODEWITHSUNDEEP
rdistancerasterterra
patch105
patch105

Reputation: 61

Distance only to Non-NA cells in Terra raster

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)

enter image description here

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)

enter image description here

This is my desired result: enter image description here

Upvotes: 2

Views: 101

Answers (3)

Andrea Manica
Andrea Manica

Reputation: 321

You could achieve what you are asking with the exclude parameter of distance if you use a single raster. To do that, you would 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)

enter image description here

However, after some benchmarking with @L Tyrone (see his excellent answer and our discussion), it turns out that a simple:

dist <- distance(r, p2)
dist <- mask(dist, r)

is computationally faster.

Upvotes: 3

L Tyrone
L Tyrone

Reputation: 7065

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)

1

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

Upvotes: 2

D.J
D.J

Reputation: 1339

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()

Upvotes: 1

Related Questions