/ sf Public

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

# Finding neighboring polygons#234

Closed
opened this issue Feb 27, 2017 · 62 comments
Closed

# Finding neighboring polygons #234

opened this issue Feb 27, 2017 · 62 comments

### tiernanmartin commented Feb 27, 2017

 I use `spdep::poly2nb` to identify neighboring polygons, but this necessitates a conversion back to `sp` - is there a way to do this without that conversion? The text was updated successfully, but these errors were encountered:

### edzer commented Feb 27, 2017

 `st_intersects` gives the sparse matrix; you'd have to filter out the diagonal elements, and set the class attribute.

### tiernanmartin commented Feb 27, 2017 • edited Loading

 Disclaimer: If Stackoverflow is the preferred venue for this sort of question please let me know and I'll post there. Below I have a reprex showing the type of polygon neighbor calculation I do regularly. The blue county has several intersecting neighbors which are identified using `sf::st_intersects`, but the southeastern-most county (light brown in color) only shares a single point with the target county and I may want to exclude it; `sp::spdep` has the `queen = FALSE` argument which allows me to easily exclude such counties. Could someone demonstrate a way to do this using `sf`-verbs? If no simple method for filtering the diagonal elements exists, could one be added? ```library(sp) library(spdep) library(tidyverse) library(sf) demo(nc, ask = FALSE, verbose = FALSE) par(mfrow = c(3, 1), mar = c(rep(0.75, 4))) # Target county plot(nc[1], col = "red", main = "Target county") plot(nc[50, "AREA"], add = TRUE) # Using sf verbs spare_mtx <- st_intersects(nc, nc) #> although coordinates are longitude/latitude, it is assumed that they are planar plot(nc[1], col = "red", main = "Neighbors: st_intersection") plot(nc[spare_mtx[[50]], "NAME"], add = TRUE) plot(nc[50, "AREA"], add = TRUE) # Using the 'rook' method from spdep::poly2nb spare_mtx_rook <- poly2nb(as(nc, "Spatial"), queen = FALSE) plot(nc[1], col = "red", main = "Neighbors: poly2nb") plot(nc[spare_mtx_rook[[50]], "NAME"], add = TRUE) plot(nc[50, "AREA"], add = TRUE)```

### edzer commented Feb 27, 2017

 Look into `st_relate`, it gives you the dimension of the intersection of the boundaries of two geometries; this is 1 for a line, 0 for a point. (It does give you a dense matrix, though; in case memory is a problem you could do `st_intersects` first, then loop through each nbh set with `st_relate`.)

### tiernanmartin commented Feb 27, 2017

 Got it. `st_relate` does return the result I needed but I would suggest that there's still room for a simpler, high-level verb to do this. Using the example above, I loop through the result of `st_intersects` which returns a list of matrices containing DE-9IM Matrix strings (which I had never seen before and take a bit of work to understand): ```spare_mtx <- st_intersects(nc, nc) spare_mtx[[50]] #> [1] 39 40 42 50 69 70 71 spare_mtx[[50]] %>% map(~nc[.x, ]) %>% map(~st_relate(x = nc[50, ], y = .x)) #> [[1]] #> [,1] #> [1,] "FF2F11212" #> #> [[2]] #> [,1] #> [1,] "FF2F11212" #> #> [[3]] #> [,1] #> [1,] "FF2F11212" #> #> [[4]] #> [,1] #> [1,] "2FFF1FFF2" #> #> [[5]] #> [,1] #> [1,] "FF2F11212" #> #> [[6]] #> [,1] #> [1,] "FF2F01212" <- Note: this is the diagonal neighbor county from the example above #> #> [[7]] #> [,1] #> [1,] "FF2F11212"``` I think a convenience function similar to `poly2nb` would be a worthwhile addition to the `sf` canon.

### tiernanmartin commented Mar 1, 2017 • edited Loading

 On a related note, I'm a bit puzzled by the following error: ```library(tidyverse) library(sf) library(purrr) demo(nc, ask = FALSE, verbose = FALSE) # Loop over nc\$geom with st_intersects() and save the results in a new list-col: nc\$INTERSECT nc %>% mutate(INTERSECT = map(.x = geom, .f = st_intersects, y = st_geometry(nc))) #> Error in mutate_impl(.data, dots): st_crs(x) == st_crs(y) is not TRUE``` In this example, it appears that the classes of the two objects are different: ```#> Browse[2]> map(list(x,y),class) #> [[1]] #> [1] "XY" "MULTIPOLYGON" "sfg" #> #> [[2]] #> [1] "sfc_MULTIPOLYGON" "sfc"``` ... and therefore `x` lacks `crs` metadata: ```#> Browse[2]> map(list(x,y),st_crs) #> [[1]] #> \$epsg #> [1] NA #> #> \$proj4string #> [1] NA #> #> attr(,"class") #> [1] "crs" #> #> [[2]] #> \$epsg #> [1] 4267 #> #> \$proj4string #> [1] "+proj=longlat +datum=NAD27 +no_defs" #> #> attr(,"class") #> [1] "crs"``` I find it confusing (and inconvenient) that a `sfc` passed to `map` results in a `sfg`, while `st_geometry(nc)` returns a `sfc`. In the example above this class mismatch causes `st_intersects` to return an error. It seems like a vignette demonstrating the recommended way of manipulating `sfc`s within a piped workflow could help resolve this type of question. Edit: For instance, this is how I solved the problem above but I have no idea if my approach uses these tools in the way they're intended to be used - some guidance would go a long way here: ```library(tidyverse) library(stringr) library(sf) library(purrr) demo(nc, ask = FALSE, verbose = FALSE) nc %>% select(NAME) %>% mutate(NB = st_intersects(geom, geom)) %>% mutate(NB_GEOMS = map(NB, ~st_geometry(nc[.x, ]))) %>% mutate(GEOM_LST = map(geom, ~st_geometry(.x) %>% st_set_crs(st_crs(nc)))) %>% mutate(RELATE = map2(.x = NB_GEOMS, .y = GEOM_LST, .f = ~st_relate(.y, .x))) %>% mutate(RELATE_LGL = map(RELATE, str_detect, pattern = "^[[:alnum:]]{4}1")) %>% mutate(NB_ROOK = map2(.x = NB, .y = RELATE_LGL, ~.x[.y])) %>% select(NAME, NB, NB_ROOK, dplyr::everything()) %>% slice(50) %>% # Rowan County is the one highlighted in the example above as_tibble() #> # A tibble: 1 Ã— 8 #> NAME NB NB_ROOK NB_GEOMS GEOM_LST #> #> 1 Rowan #> # ... with 3 more variables: RELATE , RELATE_LGL , #> # geom ```

added a commit that referenced this issue Mar 1, 2017
``` drop crs equality for geom bin ops if on arg is sfg ```
``` d9d0af3 ```
```if we call a geometry predicate like st_intersects, do not require
that x and y have identical crs when one of them (or both) is an sfg,
since sfg's don't have a crs. Adresses #234```

### edzer commented Apr 30, 2017 • edited Loading

 Your example above ```> nc %>% mutate(INTERSECT = map(.x = geom, .f = st_intersects, y = st_geometry(nc))) %>% print(n=2) Simple feature collection with 100 features and 15 fields geometry type: MULTIPOLYGON dimension: XY bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965 epsg (SRID): 4267 proj4string: +proj=longlat +datum=NAD27 +no_defs First 2 features: AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0 NWBIR74 BIR79 SID79 NWBIR79 INTERSECT geom 1 10 1364 0 19 1, 2, 18, 19 MULTIPOLYGON(((-81.47275543... 2 10 542 3 12 1, 2, 3, 18 MULTIPOLYGON(((-81.23989105...``` now works, but the simpler way to do this is `nc %>% mutate(INTERSECT = st_intersects(.))` I am not so sure I can follow your misunderstanding of `purrr::map`; although I don't think its current documentation explains well what `map` does, I believe it basically redoes `lapply`: apply a function to elements of a list, return a list with return values. Since list elements of an `sfc` are of class `sfg` (see first vignette), this is to be expected. What had you expected it to do?

### tiernanmartin commented May 2, 2017

 Thanks for adding d9d0af3 🎉 Regarding your question: I believe that my mistake was that I should have used `map2(.x = geom, .y = st_geometry(nc), .f = st_intersects)` so that both sets of elements passed to `st_intersects()` were class `sfg`. In retrospect, that example is pretty confusing and it is unclear in the post that I was using `debug()` in the second part, so thanks for even giving any thought at all.

### tiernanmartin commented May 2, 2017

 I'd like to know your thoughts (as the package developer) on the approach I demonstrated in the edit - is this the way you intend for users to step through the creation of a spatial analysis variable (such as identifying nearest neighbors)? Or do you think it's more suitable for users to compose functions that do this sort of task?

### edzer commented May 3, 2017

 With admittedly limited knowledge of `spdep`, I have never seen someone create geometries with the actual neighbours, as you seem to do in `NB_GEOMS`, rather than indices of the actual neighbours as in my example above. Why would you do so?

### tiernanmartin commented May 3, 2017

 You certainly could use the indices - it's probably more computationally efficient to do it that way. But that is not the focus of my question. My question is a bit more general: how should an `sf` user approach a multi-step spatial operation like this? Using the existing `sf` tools, finding "rook" neighbors requires at least two steps: `st_intersects()` then `st_relate()`. I demonstrated an approach that makes heavy use of `dplyr::mutate()` and `purrr::map*()`, saving intervening variables in columns. I suspect this type of operation might also be approached by writing a single function that performs all of the steps. But maybe you have a different approach altogether? I think that relatively new R-users like myself who are starting to tackle spatial problems (or spatial datasets) would benefit from a vignette demonstrating how to go about a multi-step spatial operation using `sf` tools. It looks like you've started a vignette to-do list - perhaps you might consider adding this request to your list. Thanks!

### edzer commented May 3, 2017

 I now understand it - sorry, took a while. What I think we need is to add a `pattern` variable to `st_relate`, just like `rgeos::gRelate`, which lets you select neighbours based on a pre-defined pattern, and in that case returns the indices list of `st_intersect`. You could then trivially write a function `st_rook`. Will look into.

added a commit that referenced this issue May 3, 2017
``` add parameter pattern to st_relate, address #234 ```
``` 8b61bfb ```

### edzer commented May 3, 2017

 Pls test; you should now get rook neighbors with ```st_rook = function(a, b = a) st_relate(a, b, pattern = "F***1****") nc %>% mutate(NB_ROOK = st_rook(.))```

### walkerke commented May 10, 2017

 This is fantastic! I've tested and it works. By extension, to get queen's case neighbors as well you can do: ```st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****") nc %>% mutate(NB_QUEEN = st_queen(.))``` which constitutes an improvement on using `st_intersects` for that purpose as it eliminates self-intersections. For anyone coming across this, the Wikipedia page for DE-9IM strings is quite helpful: https://en.wikipedia.org/wiki/DE-9IM.

added a commit that referenced this issue May 11, 2017
``` add queen neighbour example, #234 ```
``` 28928b4 ```

### edzer commented May 11, 2017

 Thanks @walkerke , I've added this to the documentation of `st_relate_pattern`.

### becarioprecario commented Jun 22, 2017

 Hi, I am creating a function, st2nb, to compute adjacencies from sf-objects using the the code in poly2nb() using all st_rook() and st_queen(), as explained in earlier posts in this thread. However, poly2nb() takes an additional argument called snap: snap: boundary points less than ‘snap’ distance apart are considered to indicate contiguity Is there any way to include this in the call to st_relate? Thanks! Virgilio

### edzer commented Jun 22, 2017

 I can think of two solutions: you can use `st_distance` to compute distances, and compare these with snap; this however gives you a dense matrix, set the precision, which should round away (most?) differences smaller than snap.

### becarioprecario commented Jun 22, 2017 via email

 Hi Edzer, I can think of two solutions: • you can use st_distance to compute distances, and compare these with snap; this however gives you a dense matrix, • set the precision, which should round away (most?) differences smaller than snap. Thanks for your comments. Setting the precision seems to have funny side and st_distance seems to only work point geometries. It would be doable with this but it really involves more work than using st_queen or st_rook… So I think that I will keep a warning for now. :) Best, Virgilio

### rsbivand commented Jun 22, 2017

 Somewhere I should have an object that needed non-default snap that might help. I'll try to find it next week - may we talk about this at useR!?

### becarioprecario commented Jun 22, 2017 via email

 Hi Roger, Somewhere I should have an object that needed non-default snap that might help. I'll try to find it next week - may we talk about this at useR!? Thanks! That would be very useful. I think that there is always the possibility of rewriting the existing code to deal with sf-objects instead of using st_relate, but this will require more work, I think. Best, Virgilio

### edzer commented Jun 22, 2017

 ```> p = st_polygon(list(rbind(c(0,0),c(1,1),c(0,1),c(0,0)))) > st_distance(p, p + c(2,2)) [,1] [1,] 1.414214```

### becarioprecario commented Jun 22, 2017 via email

 Hi El 22 jun 2017, a las 14:30, Edzer Pebesma ***@***.***> escribió: > p = st_polygon(list(rbind(c(0,0),c(1,1),c(0,1),c(0,0)))) > st_distance(p, p + c(2,2)) [,1] Thanks for the example. In poly2nb(), the bounding boxes of two possible neighbors are checked first to see if they overlap. If they do, then the boundaries are actually checked to see if they are neighbors or not. So, there is no need to compute all pairs of distances. But I thought that perhaps with simple features there were other (simpler) options using some of the st_XXX functions. Another option would be to add a buffer of snap/2 around each polygon and then use st_relate, right? Just thinking… Best, Virgilio

### edzer commented Jun 22, 2017

 If they do not overlap, they can be still snap distance apart... The st_buffer approach sounds good, if you insist on supporting snap.

mentioned this issue Oct 6, 2017

### Robinlovelace commented Oct 6, 2017

 I found this approach worked in a use case from sp, as documented in #210. The question is a) does this generalise to other instances where the dimension of the relation is important and b) how to make the code more user-friendly and concise? If there's no quickfire answer I can open a new issue but just wanted to check I'm not missing something obvious and easy to teach first (from an R for spatial data course where people are asking about this!): ```library(sp) g = SpatialGrid(GridTopology(c(0, 0), c(1, 1), c(3, 3))) p = as(g, "SpatialPolygons") library(sf) #> Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.2, lwgeom 2.3.3 r15473 library(tidyverse) #> Loading tidyverse: ggplot2 #> Loading tidyverse: tibble #> Loading tidyverse: tidyr #> Loading tidyverse: readr #> Loading tidyverse: purrr #> Loading tidyverse: dplyr #> Conflicts with tidy packages ---------------------------------------------- #> filter(): dplyr, stats #> lag(): dplyr, stats p_sf = st_as_sf(p) st_rook = function(a, b = a) st_relate(a, b, pattern = "F***1****") st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****") p_queen = p_sf %>% mutate(NB_ROOK = st_queen(., b = p_sf[5, ])) %>% filter(NB_ROOK == 1) nrow(p_queen) #> [1] 8 p_rook = p_sf %>% mutate(NB_ROOK = st_rook(., b = p_sf[5, ])) %>% filter(NB_ROOK == 1) nrow(p_rook) #> [1] 4 p_minDim2 = p_sf %>% filter(st_within(., y = p_sf[5, ], sparse = FALSE)) nrow(p_minDim2) #> [1] 1```

### rsbivand commented Oct 6, 2017 • edited Loading

 Notes so far on spdep's R-forge SVN site. This is very much work in progress, as `poly2nb` is pre-SF (the definition) and thus processes invalid geometries without protesting (and is much faster, even using very primitive scanning for potential neighbours). The same holds for trying st_buffer() to find distance neighbours - the vignette will have timings. Edzer has added a class (`"sgbp"` - sparse geometry binary predicate), and row.names, and what predicate made the object attributes to lists returned by `sf_relate` and friends; this will make it much easier to use sf objects inside `poly2nb`, `dnearneigh`, `knearneigh`, and the graph neighbours, as well as helper functions and methods. Because spdep (on R-forge) now suggests sf (and may import from), it has stopped checking and building on R-forge (but rgdal and rgeos build and check OK??). The report: ``````Quitting from lines 35-38 (nb_sf.Rmd) Error: processing vignette 'nb_sf.Rmd' failed with diagnostics: there is no package called 'sf' `````` But you can download the Rmd file, or checkout SVN.

### edzer commented Oct 6, 2017

 Does r-forge have GDAL >= 2.0 ?

### nuest commented Oct 7, 2017 via email

 It is possible to transfer the repository's history and match svn user names to GitHub names. I can help out if you would like to preserve these. Should the repo become part of r-spatial?

### rsbivand commented Oct 7, 2017

 Thanks for the offer of help! Is this a reasonable representation of transfer? I see that gstat is not in r-spatial, so I think I would do the same, keeping r-spatial for shared infrastructure, does that make sense?

### edzer commented Oct 7, 2017 • edited Loading

 r-spatial welcomes modern and actively maintained spatial R packages.

### nuest commented Oct 7, 2017 • edited Loading

 @rsbivand yes, that looks like a complete list of steps. This one I have used and contributed to myself a while back, includes some more steps around tags and branches, and removing large files from history (if any of these are needed).

### rsbivand commented Oct 7, 2017

 @edzer spdep isn't modern, but is actively maintained. It also invites modularisation and rationalisation, but it isn't clear that modularisation in-place (in r-spatial) is the only choice. For example, I do not use roxygen and am not convinced that it is helpful; I feel further that mechanical code coverage is not better than help page examples showing the relationship between implementation and the published books/articles with their examples. So if modern means roxygen and testthat, spdep isn't modern. If modern is without such constraints, maybe r-spatial is viable.

### edzer commented Oct 7, 2017

 With modern I was more thinking of modern, state-of-the-art methods, not of coding styles.

### rsbivand commented Oct 7, 2017

 Aha, then yes, spdep and code/packages depending on spdep are modern ... like the recent Wagner & Zeileis spatial and tree regression two-step package, and others ...

mentioned this issue Oct 10, 2017

### rsbivand commented Oct 10, 2017

 @Nowosad : may I upload a newly imported NY8 GPKG to spData? This will let me continue with the neighbour object vignette, and start a transition to importing data from spData into spdep, then spdep to github.

 @rsbivand Yes, of course. Do you prefer to make a pull request or should I give you a push access to the repository?

### rsbivand commented Oct 10, 2017

 Thanks, when I get to it, I'll make a PR.

### rsbivand commented Oct 26, 2017

 The vignette at: https://r-forge.r-project.org/scm/viewvc.php/pkg/vignettes/nb_sf.Rmd?view=markup&revision=719&root=spdep and attached now runs, see the summary section. You'll also need spData from https://github.com/rsbivand/spData; I haven't made a PR yet for this, so it isn't in Jakub Nowosad's repo. Maybe you can just try out the suggestions in the vignette summary without needing the data. I've been using the most recent sf development version. Please let me know if this is clear, and if it helps.

### edzer commented Oct 26, 2017

 When installing spdep from r-forge, after installing your spData fork, I'm seeing: ``````* creating vignettes ... ERROR Error in CPL_geos_binop(st_geometry(x), st_geometry(y), op, par, pattern, : Evaluation error: TopologyException: side location conflict at 411447.38362785219 4768469.9111850867. Quitting from lines 430-433 (nb_sf.Rmd) `````` any ideas?

### rsbivand commented Oct 26, 2017

 Checks cleanly here (two different Fedora machines). This looks like about line 269: try(NY8_sf_old_1_sgbp <- st_queen(NY8_sf_old)), the cited line numbers don't seem relevant. This has to be on the MULTIPOLYGON object, it can't be on the POINT object. Can you run it with spdep from CRAN free-standing? I'll try ... success. The error message is shown as: ``````Error in CPL_geos_binop(st_geometry(x), st_geometry(y), op, par, pattern, : Evaluation error: TopologyException: side location conflict at 411447.38362785219 4768469.9111850867. `````` in chunk 31. I've added ", silent = TRUE" to the try lines to remove the verbatim error message from the output, but R CMD build worked for me anyway when it was there.

### edzer commented Oct 26, 2017

 The error came from RANN not being installed; it now runs.

### rsbivand commented Oct 27, 2017

 Will now run without RANN too.

### rsbivand commented Nov 1, 2017

 spdep now on r-spatial on github, 0.7-2 submitted to CRAN

mentioned this issue Dec 3, 2018
mentioned this issue Sep 13, 2019
mentioned this issue Nov 19, 2019
mentioned this issue Aug 31, 2021

### rikudoukarthik commented Jan 24, 2023 • edited Loading

 This is fantastic! I've tested and it works. By extension, to get queen's case neighbors as well you can do: ```st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****") nc %>% mutate(NB_QUEEN = st_queen(.))``` which constitutes an improvement on using `st_intersects` for that purpose as it eliminates self-intersections. For anyone coming across this, the Wikipedia page for DE-9IM strings is quite helpful: https://en.wikipedia.org/wiki/DE-9IM. I created a grid using `st_make_grid()` and converted to an sf object. Now I am trying to create neighbours using the abovementioned `st_rook()` and `st_queen()` but I am getting the following error: ``````Error in xj[i, , drop = FALSE] : incorrect number of dimensions `````` I am not sure what exactly is the issue here. Am I inputting the wrong spatial object here? Or is it because the number of cells in X direction is not the same as that in Y direction? I also tried looking at the code of `st_relate()` to see if I could gather anything, but failed. I wonder if you all have some obvious solution to my issue. If not, I can share a reprex for the same. Edit: sharing reprex The necessary shapefiles can be found here ``````library(tidyverse) library(sf) country_path <- "data/in_2011/" country_file <- "India_2011" grid_sizes_km <- c(25, 50, 100, 200) grid_sizes_deg <- grid_sizes_km*1000/111111 india_sf <- st_read(dsn = country_path, layer = country_file) %>% mutate(DISTRICT = NULL) %>% st_set_crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") # creating grid g1_sf <- india_sf %>% st_make_grid(cellsize = grid_sizes_deg[1], n = (c(diff(st_bbox(india_sf)[c(1, 3)]), diff(st_bbox(india_sf)[c(2, 4)]))/grid_sizes_deg[1]) %>% ceiling()) %>% st_as_sf() %>% rename(geometry = x) # neighbours st_rook <- function(a, b = a) st_relate(a, b, pattern = "F***1****") st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****") g1_sf %>% mutate(NB_ROOK = st_rook(.), NB_QUEEN = st_queen(.)) # returns error ``````

### rsbivand commented Jan 24, 2023

 Always minimal reprex. Note that this is a very old issue, and that spdep fully supports most sf objects.

### rikudoukarthik commented Jan 24, 2023

 Always minimal reprex. Note that this is a very old issue, and that spdep fully supports most sf objects. Thanks for the very quick response. I've updated my post with a reprex. Like the OP, I would like to find a way to achieve this using only `sf`, without having to convert to and from `sp`.

### rsbivand commented Jan 24, 2023 • edited Loading

 Please avoid any tidyverse and pipes in all reprex, as they prevent access to intermediate objects, in which misunderstandings typically occur. Do not use the st_rook() etc. functions, use functions in spdep.

### rsbivand commented Jan 24, 2023

 Do not use +towgs84 ever. Do not expect spherical geometries as planar. Consider reading https://rsbivand.github.io/csds_jan23/csds_crs_workshop_230119.html if this seems hard to grasp. s2 does not behave like GEOS. See what happens in your (corrected) reprex with sf_use_s2(FALSE).

### rsbivand commented Jan 27, 2023

 @rikudoukarthik there is no problem at all in the example: ``````library(sf) country_file <- "India_2011.shp" india_sf <- st_read(country_file) st_crs(india_sf) <- "OGC:CRS84" grid_sizes_km <- c(25, 50, 100, 200) grid_sizes_deg <- grid_sizes_km*1000/111111 n <- ceiling(c(diff(st_bbox(india_sf)[c(1, 3)]), diff(st_bbox(india_sf)[c(2, 4)]))/grid_sizes_deg[1]) g1_sf <- st_make_grid(india_sf, cellsize = grid_sizes_deg[1], n = n) g1_sf st_rook <- function(a, b = a) st_relate(a, b, pattern = "F***1****") system.time(NB_ROOK <- st_rook(g1_sf)) library(spdep) system.time(NB_ROOK1 <- poly2nb(g1_sf, queen=FALSE)) all.equal(NB_ROOK1, NB_ROOK, check.attributes=FALSE) st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****") system.time(NB_QUEEN <- st_queen(g1_sf)) system.time(NB_QUEEN1 <- poly2nb(g1_sf)) all.equal(NB_QUEEN1, NB_QUEEN, check.attributes=FALSE) `````` spdep does a lot more checking (in particular for no-neighbour entities) and prepares the output object for use in spatial analysis. It does not use sp objects.

### rikudoukarthik commented Feb 19, 2023

 there is no problem at all in the example: Thanks @rsbivand for the detailed responses. I now understand the issue with using +towgs84, and have changed it as you have suggested. However, I am still facing an issue in the example. `NB_ROOK <- st_rook(g1_sf)` runs with no error message, but when I try to view it `head(NB_ROOK)` the same error I referenced in my first post pops up again (incorrect dimensions). This does not happen when using the `spdep` approach. Would you be able to shed light on what is happening here? I figure it might have to do with the fact that `st_relate()` assumes planar coordinates, which brings back the original issue. I thought `sf` now has shifted to using S2 methods in these functions, but maybe not, considering some of the discussion here. In this case, does `spdep` also use S2 methods or not? If not, although it does not return an error it can potentially cause issues, can it not?

### rsbivand commented Feb 19, 2023

 Only `c("intersects", "contains", "within", "covers", "covered_by", "disjoint", "equals", "touches")` are handled by s2, even if spherical, `c("relate_pattern", "relate")` go through GEOS, not s2. Simply stop using `st_rook()` or `st_queen()`, they have been abandoned. spdep uses s2 for spherical data.

mentioned this issue Feb 26, 2024
to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

9 participants