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

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.

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!