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
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:
- 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"))
- 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
>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
- 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).