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 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 called
- intersect polygon object
gridwith polygon of the water with
sf::st_intersectionto 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
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
Any suggestion and Idea on the workflow itself and the code is welcome.