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:
- create a grid of squared cells with desired size using
sf::st_make_grid
. this give me a grid of polygons calledgrid
. - intersect polygon object
grid
with polygon of the water withsf::st_intersection
to only pick the cells on the sea (cells partially overlapping the land are 'cutted' and not squared anymore but still present) - converting the result of the intersection (still a polygon sf object) into a raster object
- 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:
- 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?
- 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.