R cran rasterize sfc_polygon

Question

I am still trying to create a map of the Vermeille Coast in order to calculate the distance between two sampling points with the condition that the path between the two points is not crossing the land.

I am not familiar with the sf package but thanks to stack I already :

1/ bind two shapefile together (R cran: sf Sew two MULTILINESTRING/LINESTRING)

2/ created à box around to draw a polygon (Sf package: Close a polygon fom complex shape, R cran sf create a polygon from two sewed shapefiles)

3/ I am trying now to rasterize the polygon in order to apply this method (https://www.r-bloggers.com/2020/02/three-ways-to-calculate-distances-in-r/)

(data available here : https://www.dropbox.com/sh/hzsdklnmvjg4hsz/AAATHLV0pkJXDvSqyRIBlVl_a?dl=0)

To rasterize, I used the following codes:

https://cran.r-project.org/web/packages/fasterize/vignettes/using-fasterize.html https://rdrr.io/cran/fasterize/man/fasterize.html

library(sf)
library(fasterize)
library(raster)
library(dplyr)
library(devtools)
library(stars)

setwd("~/Dropbox/Data")

frenchCoast_CoteBanyuls <- st_read("coasts_subnational_France/coasts_subnational.shp")
spainCoast_CoteBanyuls <- st_read("coasts_subnational_Spain/coasts_subnational.shp")

bbox_frenchCoast_CoteBanyuls <- st_bbox(frenchCoast_CoteBanyuls) %>% 
  st_as_sfc()

polygon_frenchCoast_CoteBanyuls <- bbox_frenchCoast_CoteBanyuls %>% 
  lwgeom::st_split(frenchCoast_CoteBanyuls) %>% 
  st_collection_extract("POLYGON")

bbox_spainCoast_CoteBanyuls <- st_bbox(spainCoast_CoteBanyuls) %>% 
  st_as_sfc()

polygon_spainCoast_CoteBanyuls <- bbox_spainCoast_CoteBanyuls %>% 
  lwgeom::st_split(spainCoast_CoteBanyuls) %>% 
  st_collection_extract("POLYGON")

polyCombin <- st_union(polygon_frenchCoast_CoteBanyuls[1], polygon_spainCoast_CoteBanyuls[3])

plot(polyCombin, col = 'steelblue')

r <- raster(polyCombin, res = 1)
r <- fasterize(polyCombin, r, field = "value", fun="sum")
plot(r)

However I get the following error:

r <- raster(polyCombin, res = 1) Error in (function (classes, fdef, mtable)  :    unable to find an inherited method for function ‘raster’
for signature ‘"sfc_POLYGON"’
r <- fasterize(polyCombin, r, field = "value", fun="sum") Error in fasterize(polyCombin, r, field = "value", fun = "sum") :    Expecting
a single string value: [type=NULL; extent=0].

I then realized that in the example (https://rdrr.io/cran/fasterize/man/fasterize.html), the object to rasterize was different that mine:

library(raster)
library(fasterize)
library(sf)
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- list(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))
p3 <- list(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))
st_sf <- st_sf(value = c(1,2,3),
              geometry = st_sfc(lapply(list(p1, p2, p3), st_polygon)))

class(st_sf)
[1] "sf"         "data.frame"

in my case, I did:

class(polyCombin) [1] "sfc_POLYGON" "sfc"

So I tried to convert sfc_polygon into data.frame:

raster(as.data.frame(polyCombin), res = 1)
Error in .local(x, ...) : unused argument (res = 1)

it did not work. And I did not find any solution. I even tried :D

st_write(polyCombin, "polyCombin.shp")
polyCombin <- st_read("polyCombin.shp")
r <- raster(polyCombin, res = 1)
r <- fasterize(polyCombin, r)
plot(r)

Booh :D But it did not work.

Does someone have an idea? Thanks a lot in advance !

Charlotte

Answer 1

The problem is that you can't raster a sfc object. sfc objets are objects of sf containing ONLY the geometry. You can "add" a data.frame to an sfc object with st_sf().

Additionally, note that the CRS is important here, since your coordinates are in longitude/latitude units. If you want a raster with a specific resolution (in meters) you need to project first to a CRS defined in meters. In this case, UTM 30 zone 31N (https://spatialreference.org/ref/epsg/25831/).

I am not using your original data but an example, but the approach should still work for you:

library(giscoR)
library(sf)
library(dplyr)
library(ggplot2)


# This is mock data, use yours
polyCombin <- gisco_get_coastallines(resolution = 3, epsg = 4326) %>%
  st_geometry() %>%
  st_crop(xmin = 2.778126, ymin = 41.651276, xmax = 3.322487, ymax = 42.547747)

ggplot(polyCombin) +
  geom_sf()

# fasterize

library(raster)
library(fasterize)
# Replicate your error
# res indicate the cell size in map units
r <- raster(polyCombin, res = 100)
#> Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'raster' for signature '"sfc_POLYGON"'

class(polyCombin)
#> [1] "sfc_POLYGON" "sfc"
polyCombin
#> Geometry set for 1 feature 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 2.778126 ymin: 41.65128 xmax: 3.312999 ymax: 42.54807
#> Geodetic CRS:  WGS 84
#> POLYGON ((2.778126 41.65128, 2.781258 41.65128,...

# You need to convert to sf object (i.e. add a table)
# with a variable named var
polyCombin_df <- st_sf(var = 1, polyCombin)

class((polyCombin_df))
#> [1] "sf"         "data.frame"

polyCombin_df
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 2.778126 ymin: 41.65128 xmax: 3.312999 ymax: 42.54807
#> Geodetic CRS:  WGS 84
#>   var                     polyCombin
#> 1   1 POLYGON ((2.778126 41.65128...

# But additionally you need to project so you get the coords in meters...
# This is ETRS89 / UTM zone 31N, the right zone
st_crs(25831)$units
#> [1] "m"

polyCombin_df_t <- polyCombin_df %>% st_transform(25831)

ggplot(polyCombin_df_t) +
  geom_sf()

# And now rasterize
# res indicate the cell size in map units, in this case 100mx100m
r <- raster(polyCombin_df_t, res = 100)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
#> definition

r
#> class      : RasterLayer 
#> dimensions : 995, 443, 440785  (nrow, ncol, ncell)
#> resolution : 100, 100  (x, y)
#> extent     : 481524.9, 525824.9, 4611132, 4710632  (xmin, xmax, ymin, ymax)
#> crs        : +proj=utm +zone=31 +ellps=GRS80 +units=m +no_defs
# Ignore the warning

r <- fasterize(polyCombin_df_t, r, fun = "max")
plot(r)

^{Created on 2022-06-15 by the reprex package (v2.0.1)}

Answer 2

Try using stars and st_box to set your raster resolution:

library(stars)
library(sf)

mysf <- st_read("myshp.shp")
myras <- st_rasterize(mysf, st_as_stars(st_bbox(mysf), nx = 1000, ny = 1000, values = NA_real_))
plot(myras) #works 

Set desired raster resolution using nx and ny.