sf polygon to raster workflow to calculate least cost distance with gdistance in R

Question

I need to calclulate the shortest distance from a set of points to other sets of points in a marine area. Travelling on land must not be possible so I understand I need to calculate the least cost distance (to coherce the route to be on the water) . To do that I found the method using the R package gdistance.

Functions in gdistance work rely on raster format, and I am working with vector layers. So, in order to get to use the gdistance package i need to swithc to raster at some point.

Here's the workflow I'm trying to use so far:

1. create a grid of squared cells with desired size using sf::st_make_grid. this give me a grid of polygons called grid.
2. intersect polygon object grid with polygon of the water with sf::st_intersection to only pick the cells on the sea (cells partially overlapping the land are 'cutted' and not squared anymore but still present)
3. converting the result of the intersection (still a polygon sf object) into a raster object
4. use the converted raster in gdistance

Code wise looks like this:

library(raster)
library(sf)
# point 1
mare # polygon having only the water part of the area (by plotting it I can see land and islands shapes defined by the borders of the object mare)

>Geometry set for 1 feature 
>geometry type:  POLYGON
>dimension:      XY
>bbox:           xmin: 669911.8 ymin: 5550293 xmax: 695573.7 ymax: 5571852
>epsg (SRID):    32619
>proj4string:    +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs
>POLYGON ((682074.1 5550309, 671488.6 5553743, 6.

grid<-st_make_grid(mare, cellsize = 100, square=T) # here I specify the cell to be of 100x100 meters 

grid 

Geometry set for 55512 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 669911.8 ymin: 5550293 xmax: 695611.8 ymax: 5571893
epsg (SRID):    32619
proj4string:    +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs
First 5 geometries:
POLYGON ((669911.8 5550293, 670011.8 5550293, 6...
POLYGON ((670011.8 5550293, 670111.8 5550293, 6...
POLYGON ((670111.8 5550293, 670211.8 5550293, 6...
POLYGON ((670211.8 5550293, 670311.8 5550293, 6...
POLYGON ((670311.8 5550293, 670411.8 5550293, 6...

at this point grid is a grid of squares that cover both land and water

# point 2
grid_int<-st_intersection(grid,mare)

grid_int

Geometry set for 28037 features 
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 669911.8 ymin: 5550293 xmax: 695573.7 ymax: 5571852
epsg (SRID):    32619
proj4string:    +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs
First 5 geometries:
POLYGON ((681814.6 5550393, 681911.8 5550393, 6...
POLYGON ((681911.8 5550361, 681911.8 5550393, 6...
POLYGON ((682011.8 5550329, 682011.8 5550393, 6...
POLYGON ((682111.8 5550297, 682111.8 5550393, 6...
POLYGON ((682211.8 5550309, 682211.8 5550393, 6...


# point 3

grid_sp<-as(grid_int, "Spatial") # turn the sf object in to an sp one

r <- raster(grid_sp) # convert polygon into a raster (?)
r

class      : RasterLayer 
dimensions : 10, 10, 100  (nrow, ncol, ncell)
resolution : 2566.186, 2155.907  (x, y)
extent     : 669911.8, 695573.7, 5550293, 5571852  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

res(r)<-c(100,100)
r

class      : RasterLayer 
dimensions : 216, 257, 55512  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 669911.8, 695611.8, 5550252, 5571852  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

you can see that grid_int had 28037 features while after the conversion into a raster object I have 100 cells which have a different resolution compared to the one I created.

If I try to set the resolution of the raster to be 100x 100, I obtain a greater number of cells (i.e. 55512) than what I had in the original grid_int (i.e. 28037). My understanding is that it is just making a raster with cells over the entire bounding box, therefore I have more cells than expected, so basically I'm going back to the end of point 1 (before intersection grid had exactely 55512 features) but in a raster format. is this understanding correct?

My questions are:

1. Is this workflow correct? If so, how do I convert the 'grid_int' polygon to a raster that has the correct resolution and number of cells?
2. Is this workflow not correct? Then, I should work on a raster that cover the entire area (both land and sea) and assign different values to land and sea cells (e.g. land=0, sea=1), without doing the intersection before conversion? and then use this binary classified raster as the base to use the functions in gdistance?

Any suggestion and Idea on the workflow itself and the code is welcome.

Answer 1

This is not the workflow that I would use, but not far off. I think, in this case, that your best option for coercing "mare" to a raster is fasterize::fasterize, which incidentally is a great name for a package.

You will need a template raster to run fasterize, which you can produce using raster like you did above. Simply create a raster with the same extent as mare, then run fasterize on mare using said raster as a template.

require(fasterize)
require(raster)

template <- raster(mare, res = 100) #100m resolution
mare_raster <- fasterize(mare, template)

In this case I assume you'd like all non-water cells to be 0. This can take a bit of processing time if you're on a slower machine, but I generally find it works.

mare_raster[is.na(mare_raster)] <- 0

I am by no means an expert, please let me know if this works. Best of luck!