I'm trying to create a raster of distance to a single line feature (a 25 m contour line derived from an elevation raster). Since I want to run this in a loop and I have a ~15 million cells, I converted all cells to a shapefile of points ahead of time and saved it, so that I wouldn't need to convert the raster back to points every time I ran the loop. Only points in cells that had data (not NA) are in the shapefile (see grey areas in the picture). While all points are spaced 1 m apart, there are gaps in the point mesh due to NA's (white areas). Original raster that was converted to points spaced 1m apart. The white areas represent NA values (NoData)..

I used gDistance(), similar to this post, to find the distance between all of the points and the elevation contour line.

library(rgeos); library(raster); library(sf); library(spatial.tools)
elevation.ras = raster("elevation.tif", format = "GTiff")  
non.ras = create_blank_raster(reference_raster = elevation.ras, return_filename = F, nlayers = 1) #generate raster with footprint of elevation raster
non.ras[elevation.ras < 25] = 1 #assign 1 to all cells below 25 m
non.ras[elevation.ras >=25] = 0 #assign 0 to all cells above or at 25 m
my.line = rasterToContour(non.ras, levels = 1) #generate line at the 1/0 interface
basepoints = shapefile("basepoints.shp") #takes 15 min, 15 million points

#create distance matrix for raster
dmar = gDistance(basepoints, my.line, byid=T) #takes 60 min

This results in a matrix of 1 row and 15 million columns. Now I need to convert the matrix back to a raster, but since the area with data (non-NA) is not a perfect rectangle, I'm not sure how to get the row values to correspond with their original place in the starting raster. A straightforward approach would be to create a point mesh that includes all of the NA values, but since it takes 60 min to process just the points with data, I'm hesitant to add several million more cells.

Is there a faster alternative to gDistance()? 60 min is a long run time. I have about 25 iterations to loop through, so I would definitely like to use a faster method. I'm using Windows 10, R 3.5. My machine is decently powerful but not top-of-the-line.

  • 1
    what is create_blank_raster()? You could use sf, you'll need to st_as_sf(my.line) and use read_sf("basepoints.shp") - check out ?st_distance
    – mdsumner
    Apr 1, 2020 at 1:12
  • Oh, I forgot to include the package for the create_blank_raster() function. It's the spatial.tools package. I've added that to the post. The function generates a raster full of 0s using the footprint of a reference raster that can then be modified. I used it to simplify contour generation. read_sf speeds up the process of converting loading the shapefile into R quite a bit, but the st_distance() function takes about the same amount of time to complete. It generates nearly the same output, a single-row "units" object as opposed to a matrix of the same dimensions.
    – Mon Mo
    Apr 1, 2020 at 19:12


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.