Minimum spatial offset of points for st_intersection()
Closed this issue · 8 comments
I want to do a st_intersection() with a point (x) and polygon (y) feature.
The points are located somewhat close to each other. I already adjusted the coords a bit to have valid inputs for spatial models.
However, when doing st_intersection() I keep getting the following error.
Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) :
TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point 587176.14348371839 4784538.5250889091
Removal of problematic observations is followed by a new point pair causing this error.
Here are the coords (utm, epsg=32630) of one example pair:
| x | y |
|---|---|
| 588180.1 | 4784977 |
| 588092.5 | 4784951 |
What is the minimum required distance for a point pair to prevent this error?
Points as .shp: points.zip
Please provide a complete, reproducible example. The error message indicates a problem with your polygon, which you don't give.
library(sf)
library(magrittr)
poly <- st_read("/Users/pjs/Desktop/poly.shp")
st_read("/Users/pjs/Desktop/points/points.shp") %>%
st_intersection(poly) -> out---
Session info ---------------------------------------------------------------------------------------------------------------------------------------------------------------
setting value
version R version 3.4.0 (2017-04-21)
system x86_64, darwin16.5.0
ui RStudio (1.1.228)
language (EN)
collate en_US.UTF-8
tz Europe/Berlin
date 2017-05-17
Packages -------------------------------------------------------------------------------------------------------------------------------------------------------------------
package * version date source
DBI 0.6-1 2017-04-01 CRAN (R 3.4.0)
graphics * 3.4.0 2017-05-04 local
grDevices * 3.4.0 2017-05-04 local
grid 3.4.0 2017-05-04 local
magrittr 1.5 2014-11-22 CRAN (R 3.3.2)
methods * 3.4.0 2017-05-04 local
Rcpp 0.12.10 2017-03-19 cran (@0.12.10)
sf 0.4-3 2017-05-15 CRAN (R 3.4.0)
stats * 3.4.0 2017-05-04 local
tools 3.4.0 2017-05-04 local
udunits2 0.13 2016-11-17 CRAN (R 3.3.2)
units 0.4-4 2017-04-20 CRAN (R 3.4.0)
utils * 3.4.0 2017-05-04 local
The problem is indeed in the polygon. For me,
x = st_intersection(points, st_make_valid(poly))worked, but st_make_valid is not available in all sf installations. This also seems to work:
x = st_intersection(points, st_buffer(poly, 0))Thanks for solving that problem so quickly! Was not aware of st_make_valid() - should become handy in future problems.
I was tricked by the error message stating that the points are the problem due to "Ring Self-intersection at or near point" message.
@edzer
Whats your opinion on simplifying/rewriting such error messages to help users debugging (e.g. here to check the right feature that is causing trouble in the very first place)? The error message could also further suggest to use st_make_valid() (if that is a function which solves many invalid feature problems)?
Or do you prefer to stay with the default GEOS error message?
I'm not sure whether that would be possible.
I don't think points can ever form invalid geometries. Then, there is
> ls(2, pattern = "valid")
[1] "st_is_valid" "st_make_valid"and
> st_is_valid(poly)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE
[13] FALSE FALSE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE TRUE FALSE
[25] FALSE TRUE TRUE TRUE FALSE TRUE FALSE TRUE FALSE FALSE FALSE TRUE
[37] TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE FALSE FALSE TRUE
[49] FALSE FALSE TRUE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE TRUE
[61] TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE
[73] FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE
[85] FALSE FALSE TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE
[97] TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[109] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE
[121] TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE FALSE FALSE
[133] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[145] TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE FALSE
[157] TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE FALSE
[169] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
[181] FALSE TRUE TRUE FALSE TRUE TRUE FALSE TRUE FALSE TRUE FALSE FALSE
[193] FALSE TRUE
There were 50 or more warnings (use warnings() to see the first 50)
> warnings()
Warning messages:
1: In evalq((function (..., call. = TRUE, immediate. = FALSE, ... :
Ring Self-intersection at or near point 497592.4390760847 4797838.1264344566
2: In evalq((function (..., call. = TRUE, immediate. = FALSE, ... :
Ring Self-intersection at or near point 552757.19203570462 4724006.4269517334
# etc.and the vignette pointing out that geometries can be non-valid.
I can reproduce that your approach works on my local Ubuntu 16.04 install
> library(sf)
Linking to GEOS 3.5.1, GDAL 2.1.3, proj.4 4.9.2, lwgeom 2.3.2 r15302
It does not work on our Debian 9 Server
library(sf)
Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3, lwgeom 2.3.1 r15264
There we still get the 'Topology Exception' and a crash of the R session.
I'm trying to update to v2.3.2 and see if that works.
EDIT
I see that liblwgeom-dev is somehow not fully installed
sudo apt-get install liblwgeom-dev
Reading package lists... Done
Building dependency tree
Reading state information... Done
liblwgeom-dev is already the newest version.
0 upgraded, 0 newly installed, 0 to remove and 470 not upgraded.
1 not fully installed or removed.
Maybe this causes the crash. Investigating..
Downgraded to GEOS v3.5.1 and now it works on the Debian 9 System:
library(sf)
Linking to GEOS 3.5.1, GDAL 2.1.2, proj.4 4.9.3, lwgeom 2.3.1 r15264
I assume that the GEOS version (3.6.1, compiled from source) caused the crash. Maybe there was also some interference going on between GEOS/GDAL/lwgeom, who knows..
Can anyone recommend what to do here when one does not have access to st_make_valid(). I am getting the same issues and @pat-s. I am making use of the rnaturalearth package here for the data - apologies for not figuring out a more minimal dataset.
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
library(rnaturalearth)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
# Box to draw -------------------------------------------------------------
## Square box on projected surface
coords = list(matrix(c(-142, 61,
-138, 47,
-114, 47,
-110, 61,
-142, 61),
ncol = 2,
byrow = TRUE))
# Convert to sf -----------------------------------------------------------
outside_bc_box <- st_polygon(coords) %>%
st_sfc(crs = 4326) %>%
st_transform(3005)
# Download oceans data from rnatural earth --------------------------------
oceans <- ne_download(category = "physical", type = "ocean", scale = "medium", returnclass = "sf")
#> OGR data source with driver: ESRI Shapefile
#> Source: "C:\Users\salbers\AppData\Local\Temp\RtmpmWOtSa", layer: "ne_50m_ocean"
#> with 1 features
#> It has 3 fields
#> Integer64 fields read as strings: scalerank
oceans_prog <- oceans %>%
st_transform(3005)
# Try an st_intersection -------------------------------------------------
oceans_nw <- oceans_prog %>%
st_intersection(outside_bc_box)
#> Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)): Evaluation error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point -8826616.7648721859 -4651834.9628844205 at -8826616.7648721859 -4651834.9628844205.Created on 2018-07-11 by the reprex package (v0.2.0).
Session Info
Session info ---------------------------------------------------------------------------------------------------------
setting value
version R version 3.5.1 (2018-07-02)
system x86_64, mingw32
ui RStudio (1.2.800)
language (EN)
collate English_Canada.1252
tz America/Los_Angeles
date 2018-07-11
Packages -------------------------------------------------------------------------------------------------------------
package * version date source
assertthat 0.2.0 2017-04-11 CRAN (R 3.5.0)
backports 1.1.2 2017-12-13 CRAN (R 3.5.0)
base * 3.5.1 2018-07-02 local
bcmaps * 0.17.1 2018-03-14 CRAN (R 3.5.0)
bcmaps.rdata 0.1.5 2018-05-03 local
bindr 0.1.1 2018-03-13 CRAN (R 3.5.0)
bindrcpp * 0.2.2 2018-03-29 CRAN (R 3.5.0)
callr 2.0.4 2018-05-15 CRAN (R 3.5.0)
class 7.3-14 2015-08-30 CRAN (R 3.5.1)
classInt 0.2-3 2018-04-16 CRAN (R 3.5.0)
clipr 0.4.1 2018-06-23 CRAN (R 3.5.0)
colorspace 1.3-2 2016-12-14 CRAN (R 3.5.0)
compiler 3.5.1 2018-07-02 local
crayon 1.3.4 2017-09-16 CRAN (R 3.5.0)
datasets * 3.5.1 2018-07-02 local
DBI 1.0.0 2018-05-02 CRAN (R 3.5.0)
debugme 1.1.0 2017-10-22 CRAN (R 3.5.0)
devtools * 1.13.6 2018-06-27 CRAN (R 3.5.0)
digest 0.6.15 2018-01-28 CRAN (R 3.5.0)
dplyr * 0.7.6 2018-06-29 CRAN (R 3.5.0)
e1071 1.6-8 2017-02-02 CRAN (R 3.5.0)
evaluate 0.10.1 2017-06-24 CRAN (R 3.5.0)
ggplot2 * 3.0.0 2018-07-03 CRAN (R 3.5.0)
glue 1.2.0.9000 2018-06-19 Github (tidyverse/glue@a2c0f8b)
graphics * 3.5.1 2018-07-02 local
grDevices * 3.5.1 2018-07-02 local
grid 3.5.1 2018-07-02 local
gtable 0.2.0 2016-02-26 CRAN (R 3.5.0)
htmltools 0.3.6 2017-04-28 CRAN (R 3.5.0)
knitr 1.20 2018-02-20 CRAN (R 3.5.0)
lattice 0.20-35 2017-03-25 CRAN (R 3.5.1)
lazyeval 0.2.1 2017-10-29 CRAN (R 3.5.0)
lobstr * 0.0.0.9000 2018-04-30 Github (r-lib/lobstr@75254f1)
magrittr 1.5 2014-11-22 CRAN (R 3.5.0)
memoise 1.1.0 2017-04-21 CRAN (R 3.5.0)
methods * 3.5.1 2018-07-02 local
munsell 0.5.0 2018-06-12 CRAN (R 3.5.0)
pillar 1.2.3 2018-05-25 CRAN (R 3.5.0)
pkgconfig 2.0.1 2017-03-21 CRAN (R 3.5.0)
plyr 1.8.4 2016-06-08 CRAN (R 3.5.0)
processx 3.1.0 2018-05-15 CRAN (R 3.5.0)
purrr 0.2.5 2018-05-29 CRAN (R 3.5.0)
R6 2.2.2 2017-06-17 CRAN (R 3.5.0)
Rcpp 0.12.17 2018-05-18 CRAN (R 3.5.0)
reprex 0.2.0 2018-06-22 CRAN (R 3.5.1)
rgdal 1.3-3 2018-06-22 CRAN (R 3.5.0)
rgeos 0.3-28 2018-06-08 CRAN (R 3.5.0)
rlang 0.2.1 2018-05-30 CRAN (R 3.5.0)
rmarkdown 1.10 2018-06-11 CRAN (R 3.5.0)
rnaturalearth * 0.2.0 2018-07-11 Github (ropensci/rnaturalearth@a8fed80)
rnaturalearthdata 0.2.0 2018-07-11 local
rnaturalearthhires 0.2.0 2018-07-11 local
rprojroot 1.3-2 2018-01-03 CRAN (R 3.5.0)
rstudioapi 0.7 2017-09-07 CRAN (R 3.5.0)
scales 0.5.0.9000 2018-06-29 Github (hadley/scales@83c2e70)
sf * 0.6-3 2018-05-17 CRAN (R 3.5.0)
sp 1.3-1 2018-06-05 CRAN (R 3.5.0)
spData 0.2.9.0 2018-06-17 CRAN (R 3.5.0)
stats * 3.5.1 2018-07-02 local
stringi 1.2.3 2018-06-12 CRAN (R 3.5.0)
stringr 1.3.1 2018-05-10 CRAN (R 3.5.0)
testthat * 2.0.0 2017-12-13 CRAN (R 3.5.0)
tibble 1.4.2 2018-01-22 CRAN (R 3.5.0)
tidyselect 0.2.4 2018-02-26 CRAN (R 3.5.0)
tools 3.5.1 2018-07-02 local
units 0.6-0 2018-06-09 CRAN (R 3.5.0)
usethis * 1.3.0 2018-02-24 CRAN (R 3.5.0)
utils * 3.5.1 2018-07-02 local
whisker 0.3-2 2013-04-28 CRAN (R 3.5.0)
withr 2.1.2 2018-03-15 CRAN (R 3.5.0)