Search code examples
rrasterr-sfcran

R cran : rasterized polygon and sampling site coordinates

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:

enter image description here

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 : enter image description here

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!

Solution

  • 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)