3

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

Figure 1 - Plot showing grid cell 3 and its neighbours. Study area limit represented by dark line

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

Figure 2 - Plot showing cell 6 neighbours, and highlighting cell 1 and 15, which are not neighbour

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?

1 Answer 1

2

You probably need to set the precision of the coordinates to a non-zero value before you call st_touches. Read the documentation of st_set_precision.

2
  • Thank you Edzer, this does solve the issue when I set precision to 0.1, but when precision is set to 0.01 I get some strange results, with st_touches retrieving more boundaries than I would expect (see edit to my question). I can't exactly understand why and what st_precision changes in the object to make it work. Could you elaborate a little bit?
    – FAmorim
    Oct 28, 2019 at 13:32
  • The documentation I pointed you to points to st_as_binary, and the docs of that one has links to the info you're seeking. Oct 29, 2019 at 14:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.