distance behaviour at edges

[see24]

A while ago I created a function to calculate distance on a raster using a moving window function, it was very slow but worked well enough. Recently I realized that there is a much faster terra function for this! Before switching over I wanted to double check that the results were the same and I found that there are some differences that are largest near the edges. I wasn't sure which was correct so I plotted a band of values from each and it looks like my slow function matches the shape of the original line more closely at the edges. It is a bit hard to tell so it could be the terra function is more accurate but I thought I would share this comparison incase it is helpful.

My getDistFromSource function is included in the roads package which is on CRAN or available from GitHub here

# compare results from terra::distance with roads::getDistFromSource

library(roads)
library(magrittr)

rds <- data.frame(x = 1:1000, y = cos(1:1000/100)*500) %>% 
  sf::st_as_sf(coords = c("x", "y")) %>% 
  sf::st_union() %>% 
  sf::st_cast("LINESTRING")

rds_rast <- terra::rasterize(terra::vect(rds), 
                             terra::rast(nrows = 500, ncols = 500, 
                                         xmin = 0, xmax = 1000, 
                                         ymin = -500, ymax = 500),
                             touches = TRUE) %>% 
  terra::`crs<-`(value = "+proj=utm +zone=12")

rds_dist <- getDistFromSource(rds_rast, 400, 
                              kwidth = 1,
                              method = "terra", override = TRUE)

terra_dist <- terra::distance(rds_rast)

dif <- rds_dist - terra_dist

terra::plot(dif)

comp <- (terra_dist == 0) + (rds_dist >= 100 & rds_dist <= 102)*10 + (terra_dist >= 100 & terra_dist <= 102)*100

terra::plot(comp)

^{Created on 2024-05-16 with reprex v2.1.0}

 packageVersion("terra")
#> ‘1.7.71’

terra::gdal(lib="all")
#>    gdal     proj     geos 
#> "3.8.2"  "9.3.1" "3.11.2" 

On Windows 10

*Edited because I deprecated the getDistFromSource function but have kept the old method available for now with override = TRUE

[rhijmans]

terra::distance returns the Euclidian distance (if the crs is not lon/lat). Since getDistFromSource uses focal I assume it computes a grid-distance (you have to move from cell to cell). That would explain why it is more similar to terra::gridDist. See

g_dist = terra::gridDist(rds_rast, target=NA) 

[see24]

Ah, that makes sense! Thanks for the explanation!