sf
sf copied to clipboard
r-spatial
Loop 0 is not valid: Edge x has duplicate vertex with edge y.
This is an error that never used to arise. Is is due to an sf update or something else? Not sure how to trouble shoot this.
shapefile (strata.shp) is located: url <- "https://github.com/NOAA-EDAB/survdat/tree/master/inst/extdata/"
areaPolygon <- sf::st_read(dsn = system.file("extdata", "strata.shp",package = "survdat"), quiet = T)
sf::st_area(areaPolygon)
#> Error in s2_geography_from_wkb(x, oriented = oriented, check = check): Evaluation error: Found 3 features with invalid spherical geometry.
#> [145] Loop 0 is not valid: Edge 520 has duplicate vertex with edge 661
#> [159] Loop 0 is not valid: Edge 804 has duplicate vertex with edge 824
#> [276] Loop 0 is not valid: Edge 470 has duplicate vertex with edge 482.
Created on 2021-08-17 by the reprex package (v2.0.1)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
Created on 2021-08-17 by the reprex package (v2.0.1)
Session Info
-
Session info -------------------------------------------------------------------------- setting value
version R version 4.0.5 (2021-03-31) os Windows 10 x64
system x86_64, mingw32
ui RStudio
language (EN)
collate English_United States.1252
ctype English_United States.1252
tz America/New_York
date 2021-08-17 -
Packages ------------------------------------------------------------------------------ package * version date lib source
assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.3) class 7.3-19 2021-05-03 [1] CRAN (R 4.0.5) classInt 0.4-3 2020-04-07 [1] CRAN (R 4.0.3) cli 3.0.1 2021-07-17 [1] CRAN (R 4.0.5) clipr 0.7.1 2020-10-08 [1] CRAN (R 4.0.3) crayon 1.4.1 2021-02-08 [1] CRAN (R 4.0.4) data.table 1.14.0 2021-02-21 [1] CRAN (R 4.0.4) DBI 1.1.1 2021-01-15 [1] CRAN (R 4.0.3) dplyr 1.0.7 2021-06-18 [1] CRAN (R 4.0.5) e1071 1.7-8 2021-07-28 [1] CRAN (R 4.0.5) ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.0.5) fansi 0.5.0 2021-05-25 [1] CRAN (R 4.0.5) generics 0.1.0 2020-10-31 [1] CRAN (R 4.0.3) glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.3) KernSmooth 2.23-20 2021-05-03 [1] CRAN (R 4.0.5) lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.4) magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.0.3) pillar 1.6.2 2021-07-29 [1] CRAN (R 4.0.5) pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.3) proxy 0.4-26 2021-06-07 [1] CRAN (R 4.0.5) purrr 0.3.4 2020-04-17 [1] CRAN (R 4.0.3) R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.3) Rcpp 1.0.7 2021-07-07 [1] CRAN (R 4.0.5) rlang 0.4.10 2020-12-30 [1] CRAN (R 4.0.3) rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.0.3) s2 1.0.6 2021-06-17 [1] CRAN (R 4.0.5) sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.3) sf 1.0-2 2021-07-26 [1] CRAN (R 4.0.5) survdat 1.0 2021-07-16 [1] local
tibble 3.1.3 2021-07-23 [1] CRAN (R 4.0.5) tidyselect 1.1.1 2021-04-30 [1] CRAN (R 4.0.5) units 0.7-2 2021-06-08 [1] CRAN (R 4.0.5) utf8 1.2.2 2021-07-24 [1] CRAN (R 4.0.5) vctrs 0.3.8 2021-04-29 [1] CRAN (R 4.0.5) withr 2.4.2 2021-04-18 [1] CRAN (R 4.0.5) wk 0.5.0 2021-07-13 [1] CRAN (R 4.0.5)
[1] C:/~/R/win-library/4.0 [2] C:/Program Files/R/R-4.0.5/library
It would be great if you could provide the data set that generated this, if needed off-line (email), to help us understand what causes this.
The solution to this is to first set
sf_use_s2(FALSE)
Also appears with rworldmap
library(rworldmap)
data(countryExData)
p <- st_sf(geometry= st_sfc(st_polygon(list(matrix(c(
3, 40,
3, 51,
25, 51,
25, 40,
3, 40),5,2, byrow=T)))), crs="+proj=latlong +datum=WGS84")
cntr <- st_intersection(st_as_sf(countriesLow), p)
But the workaround helps.
In case it's useful there's another example posted at http://www.math.mcmaster.ca/bolker/misc/wwf_terr_ecos.shp
library(spdep)
shp <- as(sf::st_read("wwf_terr_ecos.shp"), "Spatial")
## sf_use_s2(FALSE)
p <- poly2nb(shp)
If this is a problem with the data set, could you point me to ways of diagnosing/fixing it ...?
@bbolker please contribute the whole set of *.shp, *.shx, *.dbf at least, probably also *.prj. Four files are needed if the position data have a CRS. Note that spdep::poly2nb() has accepted sf objects for years. The shapefile is large, probably old, probably therefore poorly specified, and subject to all the manifold drawbacks of a deprecated format.
Will do, thanks. I got it second-hand, but getting the latest (??) version from here: https://www.sciencebase.gov/catalog/item/508fece8e4b0a1b43c29ca22 (and the version from https://www.worldwildlife.org/publications/terrestrial-ecoregions-of-the-world) still throws the same error.
shp <- sf::st_read("wwf2/wwf_terr_ecos.shp")
p <- poly2nb(shp)
Thanks! As an alternative you might try
p <- poly2nb(st_make_valid(shp))
if you want to stay on a round rather than a flat Earth. We explain here why many polygons cannot be valid on S2 and R2 simultaneously.
> shp <- sf::st_read("official/wwf_terr_ecos.shp")
Reading layer `wwf_terr_ecos' from data source
`/home/rsb/tmp/bigshape/official/wwf_terr_ecos.shp' using driver `ESRI Shapefile'
Simple feature collection with 14458 features and 21 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -180 ymin: -89.89197 xmax: 180 ymax: 83.62313
Geodetic CRS: WGS 84
> p <- spdep::poly2nb(shp)
Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, :
Loop 15 is not valid: Edge 0 has duplicate vertex with edge 4
> table(sf::st_is_valid(shp))
FALSE TRUE
1 14457
> shp[!sf::st_is_valid(shp),]
Simple feature collection with 1 feature and 21 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -63.00804 ymin: 1.112058 xmax: -50.96594 ymax: 8.663469
Geodetic CRS: WGS 84
OBJECTID AREA PERIMETER ECO_NAME REALM BIOME ECO_NUM
1526 1501 475923.6 75.484 Guianan moist forests NT 1 25
ECO_ID ECO_SYM GBL_STAT G200_REGIO G200_NUM G200_BIOME
1526 60125 70 3 Guianan Moist Forests 42 1
G200_STAT Shape_Leng Shape_Area area_km2 eco_code PER_area PER_area_1
1526 0 75.31988 38.78856 476136 NT0125 0 0
PER_area_2 geometry
1526 0 POLYGON ((-62.54998 7.93176...
> p <- spdep::poly2nb(sf::st_make_valid(shp))
> p
Neighbour list object:
Number of regions: 14458
Number of nonzero links: 16076
Percentage nonzero weights: 0.007690624
Average number of links: 1.11191
8601 regions with no links:
5 8 9 18 20 23 24 25 26 27 28 29 31 32 33 34 35 36 37 39 44 45 46 47 48
50 51 52 53 58 59 61 62 64 65 66 69 70 72 74 76 77 79 83 87 89 90 91 92
93 94 95 96 99 100 101 102 104 105 106 107 108 110 113 114 115 116 125
150 151 153 184 187 188 189 190 191 201 204 206 208 232 250 254 255 264
268 275 279 296 302 303 306 311 317 319 320 324 325 326 335 339 342 343
...
I assume that spdep::poly2nb() isn't the intended use. The polygons do not seem to be a tessellation. More than half have no neighbour for many values of snap=, possibly islands of course, but a bit meaningless, and some overlap.
Thanks, that's very useful! I'm feeling my way around a bit with this (cargo-culting), will talk to my more spatially savvy collaborator. You've been immensely helpful.
Dear Prof. @rsbivand @edzer , It seems that both "sf_use_s2(F)", and "st_make_valid(shp)" do not work for me. How can I fix the following errors?
Best wishes, Shunsen
county_shape2 <- read_sf("ssss.shp") nb2 <- poly2nb(county_shape2) Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 127 has duplicate vertex with edge 129 weight_matrix2 <- nb2listw(nb2, style = "W") Error in nb2listw(nb2, style = "W") : Empty neighbour sets found
st_is_longlat(county_shape2) [1] TRUE sf_use_s2() [1] TRUE sf_use_s2(F) Spherical geometry (s2) switched off nb2 <- poly2nb(county_shape2) although coordinates are longitude/latitude, st_intersects assumes that they are planar weight_matrix2 <- nb2listw(nb2, style = "W") Error in nb2listw(nb2, style = "W") : Empty neighbour sets found nb2 <- poly2nb(st_make_valid(county_shape2)) although coordinates are longitude/latitude, st_intersects assumes that they are planar There were 18 warnings (use warnings() to see them) warnings() warnings： 1: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 2: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 3: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 4: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 5: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 6: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 7: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 8: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 9: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 10: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 11: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 12: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 13: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 14: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 15: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 16: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 17: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained 18: In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) : only first part of geometrycollection is retained weight_matrix2 <- nb2listw(nb2, style = "W") Error in nb2listw(nb2, style = "W") : Empty neighbour sets found
nb2 <- poly2nb(st_make_valid(county_shape2)) although coordinates are longitude/latitude, st_intersects assumes that they are planar There were 18 warnings (use warnings() to see them) weight_matrix2 <- nb2listw(nb2, style = "W") Error in nb2listw(nb2, style = "W") : Empty neighbour sets found sf_use_s2(T) Spherical geometry (s2) switched on nb2 <- poly2nb(st_make_valid(county_shape2)) Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 0 crosses edge 46 weight_matrix2 <- nb2listw(nb2, style = "W") Error in nb2listw(nb2, style = "W") : Empty neighbour sets found
Reproducible example needed. Did you examine the output nb object to check for no-neighbour observations? Did you try a different snap= value? This is not an sf issue.
Dear Prof. @rsbivand , Thank your feedback.
I have tried the snap=value, and check for no-neighbour observations. here are the results, and I attached the reporducible example and the .shp file. So if it is not the sf issue, could it be the issue of the .shp file? if it is, how can I fix it? Because only using st_make_valid did not work well.
Best wishes,
results：
county_shape2 <- read_sf("中华人民共和国.shp") #county_shape2 <- county_shape2[!is.na(county_shape2$name) ,] sf_use_s2() [1] TRUE eps <- sqrt(2.220446e-22) table(st_is_valid(county_shape2))
FALSE TRUE 20 15
county_shape_valid <- st_make_valid(county_shape2) table(st_is_valid(county_shape_valid))
FALSE TRUE 5 30
summary(county_shape_valid) adcode name center centroid
Min. :110000 Length:35 Min. : NA Min. : NA
1st Qu.:312500 Class :character 1st Qu.: NA 1st Qu.: NA
Median :425000 Mode :character Median : NA Median : NA
Mean :423235 Mean :NaN Mean :NaN
3rd Qu.:537500 3rd Qu.: NA 3rd Qu.: NA
Max. :820000 Max. : NA Max. : NA
NA's :1 NA's :35 NA's :35
childrenNu level parent subFeature
Min. : 0.00 Length:35 Min. : NA Min. : 0.00
1st Qu.:10.25 Class :character 1st Qu.: NA 1st Qu.: 8.25
Median :14.00 Mode :character Median : NA Median :16.50
Mean :13.97 Mean :NaN Mean :16.50
3rd Qu.:16.00 3rd Qu.: NA 3rd Qu.:24.75
Max. :38.00 Max. : NA Max. :33.00
NA's :1 NA's :35 NA's :1
acroutes adchar geometry Min. : NA Length:35 MULTIPOLYGON :35
1st Qu.: NA Class :character epsg:4326 : 0
Median : NA Mode :character +proj=long...: 0
Mean :NaN
3rd Qu.: NA
Max. : NA
NA's :35
nb2 <- poly2nb(county_shape2,snap = eps) Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 127 has duplicate vertex with edge 129 weight_matrix2 <- nb2listw(nb2, style = "W")
Reproducible example needed. Did you examine the output nb object to check for no-neighbour observations? Did you try a different snap= value? This is not an sf issue.
You only included a *.shp file. A shapefile is at least 3 files, *.shx, *.dbf, and preferably *.prj, and *.cpg if coming from a Windows codepage. For me, only UTF-8 is legible, so please change the basename of the files to ASCII. The single *.shp file can be supplemented by a guessed *.shx file by SHAPE_RESTORE_SHX=YES ogrinfo x.shp. On that basis:
library(sf)
x <- read_sf("x.shp")
st_crs(x) <- "EPSG:4326"
sf_use_s2(FALSE) # the input geometries are bad, and wrong spherically
x1 <- st_make_valid(x)
all(st_is_valid(x1))
library(spdep)
nb <- poly2nb(x1) # snap not required
nb
plot(st_geometry(x1), border="grey")
plot(nb, st_geometry(x1), add=TRUE)
where the three entities without neighbours are islands, so do not share boundaries with other entities. Consider reporting back to the provider of the original shapefile that the polygons are not proper spherical geometries. See also https://r-spatial.org/book/04-Spherical.html.