R: Should I transform my Raster to spatial objects here?

Question

This question is a follow-up based on some of the advice provided in the answers and comments at Converting raster to spatial points without loss of information in R.

First, to provide a bit of context, I have a geographic area of interest (NSW in Australia), and I also have a shapefile that contains the locations of drainages in the same area. I have a raster which contains pixels in the area that I am interested in (this raster is stored in a variable called ndvi_median).

My goal is to create a raster where the pixels are the area of interest, and with the raster attribute being the shortest distance (in m, although km is also fine) between each pixel and the nearest drainage. Fortunately, there was a similar question asked on this stack exchange before: Finding distance between raster pixels and line features in R.

The only difference is that my shapefiles & rasters are in lat-long, but I want the distances in metres. I have read from another stack exchange post that we should project to another CRS - I'll be using +proj=lcc +lat_1=-30.75 +lat_2=-35.75 +lat_0=-33.25 +lon_0=147 +x_0=9300000 +y_0=4500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs for this task. So, my general approach is:

1. Reproject where necessary using new CRS:

> ndvi_median # contains pixels in the geographic area of interest
class      : RasterLayer 
dimensions : 1122, 181, 203082  (nrow, ncol, ncell)
resolution : 0.09978817, 0.008333333  (x, y)
extent     : 140.9993, 159.0609, -37.50702, -28.15702  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=GRS80 +no_defs 
source     : C:/Users/thoma/AppData/Local/Temp/RtmpWGYBpT/raster/r_tmp_2021-05-18_120505_15536_85040.grd 
names      : layer 
values     : -1966, 9009.5  (min, max)

> GeoDat_NSW ## Contains information on drainages
class       : SpatialLinesDataFrame 
features    : 1868 
extent      : 140.9993, 153.5635, -37.50702, -28.15702  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
variables   : 12
names       : FNODE_, TNODE_, LPOLY_, RPOLY_,            LENGTH, AUS25DGD_,   ID,    FEAT_CODE,              NAME, PERENNIAL,   Q_INFO,        UFI 
min values  :   5952,   5952,      0,      0, 0.003601388621079,      6150, 6150,        canal, ABERCROMBIE RIVER,         0, BI000001, DJ00010416 
max values  :   9039,   8972,      0,      0,  3.84625259743069,      8737, 8737, watercours_l,      YOWRIE RIVER,         2, BJ000004, DJ00013003 

>Spat_points_NSW<-projectRaster(ndvi_median,crs="+proj=lcc +lat_1=-30.75 +lat_2=-35.75 +lat_0=-33.25 +lon_0=147 +x_0=9300000 +y_0=4500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
>Spat_points_NSW_2<-as(Spat_points_NSW,"SpatialPoints")
>GeoDat_Line<-as(GeoDat_NSW,"SpatialLines")
>GeoDat_Line_transformed<-spTransform(GeoDat_Line,CRS("+proj=lcc +lat_1=-30.75 +lat_2=-35.75 +lat_0=-33.25 +lon_0=147 +x_0=9300000 +y_0=4500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))

2. Now that everything has been reprojected using the new CRS, I now apply the method suggested at Finding distance between raster pixels and line features in R to find the distance between cells and the nearest drainage line. In fact, this is a key reason as to why I felt the need to convert everything to sp objects (which are required in the arguments of the gdistance function).

>Dist_drainage<-gDistance(Spat_points_NSW_2,GeoDat_Line_transformed,byid=TRUE) # finds distance between each cell and every drainage line 
>Min_Dist<-apply(Dist_drainage,2,min) # finds the minimum distance for each cell
>Spat_points_NSW_3<-Spat_points_NSW
>Spat_points_NSW_3[]=Min_Dist # Effectively creating a raster with the distance information as the raster attributes

3. Now, I just need to convert the coordinates in Spat_points_NSW_3 back into lat-long.

# Reproject to lat-long
>Spat_points_NSW_4<-projectRaster(Spat_points_NSW_3,crs=crs(ndvi_median))

Unfortunately, this method doesn't completely work because in the line below, I get an error

> Min_Dist<-apply(Dist_drainage,2,min)
Error: cannot allocate vector of size 1.7 Gb

This brings me to ask whether there is an alternative way to approach this problem? I saw a promising solution using geosphere::dist2Line at Calculate distance between points and nearest polygon in R, but unfortunately the method takes too long to run (although I only stopped it after running it for roughly 40mins - I can try to leave it running overnight if it seems to be the most optimal approach).

Answer 1

It would be very helpful if you provided a minimal, self-contained, reproducible example.

Here is an approach with terra that works with the example data. I have not tried it with larger data sets. But it rasterizes the cells touched by the line, and then computes the distance to the center of these cells; so it introduces some error.

library(terra)
r <- rast(ncol=36,nrow=18, crs="+proj=longlat +datum=WGS84")
v <- vect(matrix(c(-120, 47, -137, 9, -124, -12, -129, -41), ncol=2, byrow=TRUE), type="lines", crs=crs(r))

d <- distance(r, v) 
plot(d / 1000)
lines(v)

Contrary to popular belief, transforming lon/lat data to a planar CRS is not a good approach to measure distance (or area), unless you are using a planar CRS that does not distort distance (or is equal area). But it is true that geos::gDistance requires it. I will implement the geosphere method in terra to make it run (hopefully much) faster.