I am trying to get the total length of all the shared boundaries between cells in a grid (following this post). However, I noticed that some of the objects have no associated boundaries, even though they do exist.
Because the issue might be related to my specific objects, you can download the initial object (study.area) here. The code bellow allows you to create the grid and the boundaries between cells.
require(sf)
require(ggplot2)
load("file_path/study_area.RData") #load study area into R, repace "file_path" by your directory
# create hexagonal grid within study_area
cell_diameter <- sqrt(2 * 10000 / sqrt(3)) # set grid area to 10000m (1ha)
grid <- st_make_grid(study.area, cellsize = cell_diameter, square = FALSE) # create grid using "st_make grid" in sf package
grid <- st_sf(grid) # convert to sf
grid$id <- 1:length(grid$grid) # add na indentifying id to each grid cell
Touching_List <- st_touches(grid) # create list identifying all the neighbouring grids
By calling Touching_List you can see that, for example, list object 1, 2 and 3 are empty, that said, they have no boundaries.
Touching_List
Sparse geometry binary predicate list of length 16480, where the predicate was `touches'
first 10 elements:
1: (empty)
2: (empty)
3: (empty)
4: 6
5: 7
6: 4
7: 5
8: 9
9: 8
10: (empty)
So, I decided to have a better look on this by plotting grid cell 3 and its neighbours.
# build demonstration plot
ggplot(data = grid[c(3, 1, 6, 4),]) +
geom_sf(fill = c("orangered", "peachpuff1", "peachpuff1", "peachpuff1")) +
geom_sf(data = shp.avz, fill = NA, size = 1, colour = "black") +
geom_sf_label(aes(label = id)) +
coord_sf(xlim = c(-27570, -27325), ylim = c(6935, 7260), expand = FALSE) +
theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
As you can see in the picture, cell 3 does have neighbours although this is not shown in Touching_List. Finally, the last step is to check whether grid cell 3 actually appears as a neighbour on cell 6, 4 and 1.
Touching_List[c(1,4,6)]
[[1]]
integer(0)
[[2]]
[1] 6
[[3]]
[1] 4
And it does not. Finally, I noticed that this only happens with cells located at the limit of my study area (see Figure 1).
EDIT 1
This Edit uses st_precision as suggested in the answer provided by Edzer Pebesma. Although setting precision to 0.1 does seem to retrieve the expected results, setting it for instance to 0.01 returns unexpected boundaries.
# Precision set to 0.1
st_precision(grid) <- 0.1
Touching_List <- st_touches(grid)
head(Touching_List)
[[1]]
[1] 2 3 4
[[2]]
[1] 1 4 5
[[3]]
[1] 1 4 6
[[4]]
[1] 1 2 3 5 6 7
[[5]]
[1] 2 4 7 8
[[6]]
[1] 3 4 7 10
As you can see in Figure 1, st_touches now returns the expected neighbours for cell 3 (cells 1, 4 and 6).
Now, let's instead set precision to 0.01.
# Precision set to 0.01
st_precision(grid) <- 0.01
Touching_List <- st_touches(grid)
head(Touching_List)
[[1]]
[1] 2 3 4 6
[[2]]
[1] 1 4 5 7
[[3]]
[1] 1 4 6
[[4]]
[1] 1 2 3 5 6 7 10
[[5]]
[1] 2 4 7 8 11
[[6]]
[1] 1 3 4 7 10 15
In this case, we get more neighbours than expected. Let us take cell 6 as an example, as you see in the previous list st_touches returns cells 1, 3, 4, 7, 10 and 15 as neighbours. Plotting this cells
ggplot(data = grid[c(6, 1, 3, 4, 7, 10, 15),]) +
geom_sf(fill = c("orangered", "yellow", "peachpuff1", "peachpuff1", "peachpuff1", "peachpuff1", "yellow")) +
geom_sf(data = shp.avz, fill = NA, size = 1, colour = "black") +
# geom_sf(data = shp.avz, aes()) +
geom_sf_label(aes(label = id)) +
coord_sf(xlim = c(-27570, -27200), ylim = c(6935, 7450), expand = FALSE) +
theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
From the picture, cell 1 and 15 are clearly not neighbours of cell 6, although they are identified by st_touches when precision is set to 0.01. Any thoughts on why this happens?
