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). .
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.
?st_distance
read_sf
speeds up the process of converting loading the shapefile into R quite a bit, but thest_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.