R cran : rasterized polygon and sampling site coordinates

Question

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

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

2/ We created a box around to draw a polygon (Sf package: Close a polygon fom complex shape)

3/ We did rasterize the polygon (R cran rasterize sfc_polygon)

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

library(sf)
library(fasterize)
library(raster)
library(dplyr)
library(tidyverse)

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

spainurl <- "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:coasts_subnational&outputformat=SHAPE-ZIP&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Emrgid_1%3C%2FPropertyName%3E%3CLiteral%3E3417%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E"
download.file(spainurl, "spain.zip", mode = "wb")
unzip("spain.zip", exdir = "spain", junkpaths = TRUE)

franceurl <- "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:coasts_subnational&outputformat=SHAPE-ZIP&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Emrgid_1%3C%2FPropertyName%3E%3CLiteral%3E19888%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E"
download.file(franceurl, "france.zip", mode = "wb")
unzip("france.zip", exdir = "france", junkpaths = TRUE)

spainCoast_CoteBanyuls <- list.files("spain",
                                     pattern = "shp$",
                                     full.names = TRUE) %>% st_read()
frenchCoast_CoteBanyuls <- list.files("france",
                                      pattern = "shp$",
                                      full.names = TRUE) %>% st_read()


lines_spain <- st_geometry(spainCoast_CoteBanyuls) %>% st_cast("LINESTRING")
spainCoast_l <- st_sf(n = as.character(seq_len(length(lines_spain))), lines_spain)

lines_france <- st_geometry(frenchCoast_CoteBanyuls) %>% st_cast("LINESTRING")
franceCoast_l <- st_sf(n = as.character(seq_len(length(lines_france))), lines_france)

spainmax <- spainCoast_l[which.max(st_length(spainCoast_l)), ]
spainrest <- spainCoast_l[-which.max(st_length(spainCoast_l)), ]

joined <- c(st_geometry(spainmax), st_geometry(franceCoast_l)) %>%
  st_union()

join_end <- st_union(joined, st_geometry(spainrest))

bbox_all <- st_bbox(joined) %>%
  st_as_sfc()

polygon_joined <- bbox_all %>%
  lwgeom::st_split(join_end) %>%
  st_collection_extract("POLYGON")

#Polygons on position 2 and 3 need to be removed (visual inspection)
polygon_end <- polygon_joined[2] # define land as polygone and not sea

polyCombin_df <- st_sf(var = 1, polygon_end)
class(polyCombin_df)
st_crs(25831)$units
polyCombin_df_t <- polyCombin_df %>% st_transform(25831)

We got that:

and we did rasterize the polygon :

r <- raster(polyCombin_df_t, res = 100)
r <- fasterize(polyCombin_df_t, r, fun = "max")
par(mar=c(1,1,1,1))
plot(r)

So far it seems to work :

4/ We want now to add 3 sampling sites coordinates along the coast using the function points, in order to apply the following method to calculate the distance between sampling points : (https://www.r-bloggers.com/2020/02/three-ways-to-calculate-distances-in-r/)

# sites
site_random <- matrix(data = c(3.164887 , 3.123969 , 3.158125 , 3.160378, 42.402158, 
                42.521957, 42.475956, 42.461188), ncol = 2)

site_random <- data.frame(site_random) 
points(site_random$X1, site_random$X2, pch = 19)

However, this does not work and no sample site is displayed. Is this due to the scale of the graph?

Thank you in advance for your help!

Try to convert `site_random` to an `sf` object setting the crs to 4326 and then transform it to the same crs than your raster 25830. Then you should use `plot(st_geometry(points_sf), add = TRUE)` — dieghernan, Jun 22 '22 at 16:23
The data in not (anymore?) available. However you should maybe try terra (https://www.rspatial.org/terra/spatial/2-spatialdata.html#simple-representation-of-spatial-data). — defuneste, Jun 22 '22 at 16:24
Core problem here IMHO is that the sampling points are in longitude latitude coordinates whereas the raster is not (note the values of the axis in the plot). So the coordinates of the raster and the coordinates of the points needs to be in the same CRS — dieghernan, Jun 22 '22 at 16:29

score 1 · Accepted Answer · answered Jun 23 '22 at 07:52

See how to aling the coordinates of your points to the CRS of the raster:

library(sf)
library(fasterize)
library(raster)
library(dplyr)
library(tidyverse)

spainurl <- "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:coasts_subnational&outputformat=SHAPE-ZIP&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Emrgid_1%3C%2FPropertyName%3E%3CLiteral%3E3417%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E"
download.file(spainurl, "spain.zip", mode = "wb")
unzip("spain.zip", exdir = "spain", junkpaths = TRUE)

franceurl <- "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:coasts_subnational&outputformat=SHAPE-ZIP&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Emrgid_1%3C%2FPropertyName%3E%3CLiteral%3E19888%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E"
download.file(franceurl, "france.zip", mode = "wb")
unzip("france.zip", exdir = "france", junkpaths = TRUE)

spainCoast_CoteBanyuls <- list.files("spain",
  pattern = "shp$",
  full.names = TRUE
) %>% st_read()

frenchCoast_CoteBanyuls <- list.files("france",
  pattern = "shp$",
  full.names = TRUE
) %>% st_read()
lines_spain <- st_geometry(spainCoast_CoteBanyuls) %>% st_cast("LINESTRING")
spainCoast_l <- st_sf(n = as.character(seq_len(length(lines_spain))), lines_spain)

lines_france <- st_geometry(frenchCoast_CoteBanyuls) %>% st_cast("LINESTRING")
franceCoast_l <- st_sf(n = as.character(seq_len(length(lines_france))), lines_france)

spainmax <- spainCoast_l[which.max(st_length(spainCoast_l)), ]
spainrest <- spainCoast_l[-which.max(st_length(spainCoast_l)), ]

joined <- c(st_geometry(spainmax), st_geometry(franceCoast_l)) %>%
  st_union()

join_end <- st_union(joined, st_geometry(spainrest))

bbox_all <- st_bbox(joined) %>%
  st_as_sfc()

polygon_joined <- bbox_all %>%
  lwgeom::st_split(join_end) %>%
  st_collection_extract("POLYGON")

# Polygons on position 2 and 3 need to be removed (visual inspection)
polygon_end <- polygon_joined[-c(2, 3)] # define land as polygone and not sea

polyCombin_df <- st_sf(var = 1, polygon_end)
class(polyCombin_df)
#> [1] "sf"         "data.frame"
st_crs(25831)$units
#> [1] "m"
polyCombin_df_t <- polyCombin_df %>% st_transform(25831)


ggplot(polyCombin_df_t) +
  geom_sf()


r <- raster(polyCombin_df_t, res = 100)
r <- fasterize(polyCombin_df_t, r, fun = "max")



# sites
site_random <- matrix(data = c(
  3.164887, 3.123969, 3.158125, 3.160378, 42.402158,
  42.521957, 42.475956, 42.461188
), ncol = 2)

site_random <- data.frame(site_random)

site_random_sf <- st_as_sf(site_random, coords = c("X1", "X2"), crs = 4326) %>%
  st_transform(25831)


plot(r)
plot(st_geometry(site_random_sf), add = TRUE)

# If you need as matrix points
site_random_t <- st_coordinates(site_random_sf)[, 1:2] %>% as.data.frame()
names(site_random_t) <- names(site_random)

site_random_t
#>         X1      X2
#> 1 513569.2 4694442
#> 2 510182.5 4707739
#> 3 512997.5 4702635
#> 4 513185.8 4700996

plot(r)
points(site_random_t$X1, site_random_t$X2, pch = 19)

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

Thank you very much for your answer to my problem. However, when I export the plot and look at it precisely, it seems that the points of the sites are not placed with great precision along the coast. It seems to me that this is due to the rasterization of the coastline. Is it possible to increase the pixelization while rasterizing? — Florian B., Jun 29 '22 at 11:12
Yes, use for example `raster(polyCombin_df_t, res = 50)` `res` parameter controls the size of the raster pixels — dieghernan, Jun 30 '22 at 07:23

R cran : rasterized polygon and sampling site coordinates

1 Answers1