Closed bcaradima closed 5 years ago
edzer
commented
6 years ago This is what happens when doing this kind of calculations with floating point doubles. Setting precision before doing the computations might help; read the help of st_set_precision and st_as_binary.
bcaradima
commented
6 years ago Thank you, Edzer. If I understand the help documents correctly, lowering the precision to an acceptable value (considering the spatial scale of the analysis) would allow st_intersect (and other geometric operations) to round the input coordinates of a point and line 3.501388e-11m apart and to confirm that they do intersect. Using the same data, setting the precision to 1e-05m yields the following:
# see for explanation: https://r-spatial.github.io/sf/reference/st_as_binary.html
x <- st_distance(pt.nearest, reach.nearest)
precision <- 1e-05
round(x*precision)/precision
pt.nearest <- st_set_precision(pt.nearest, precision = 1e-05)
reach.nearest <- st_set_precision(reach.nearest, precision = 1e-05)
st_distance(pt.nearest, reach.nearest)
st_intersects(pt.nearest, reach.nearest, sparse=FALSE)Would it be a worthwhile feature set the precision when using st_read?
edzer
commented
6 years ago I don't think you want that much rounding:
> st_set_precision(pt.nearest, precision = 1e-05) %>% st_as_binary %>% st_as_sfc
Geometry set for 1 feature
geometry type: POINT
dimension: XY
bbox: xmin: 6e+05 ymin: 2e+05 xmax: 6e+05 ymax: 2e+05
epsg (SRID): NA
proj4string: NA
POINT (6e+05 2e+05)My mistake! My coordinates are in meters, so if the nearest site and nearest line coordinates differ at three or four decimal places then that should be an acceptable intersection. So I should set a scale argument of 1000 (to specify 3 decimal places of precision). If the scale argument equals zero, the distance between the two points is 3.501388e-11 meters. But if I use a scale of 1000, the distance increases:
> p <- 1000
> pt.nearest <- st_set_precision(pt.nearest, precision = p)
> reach.nearest <- st_set_precision(reach.nearest, precision = p)
> st_distance(pt.nearest, reach.nearest)
Units: m
[,1]
[1,] 0.0001278386
> st_intersects(pt.nearest, reach.nearest, sparse=FALSE)
[,1]
[1,] FALSE
>
> st_set_precision(pt.nearest, precision = p) %>% st_as_binary %>% st_as_sfc
Geometry set for 1 feature
geometry type: POINT
dimension: XY
bbox: xmin: 605979 ymin: 220842.5 xmax: 605979 ymax: 220842.5
epsg (SRID): NA
proj4string: NA
POINT (605979 220842.5)That doesn't really make sense, unless I'm missing something.
edzer
commented
6 years ago Maybe you should set precisions before computing pt.nearest etc.
bcaradima
commented
6 years ago Setting the scale arguments of the example data gives an intersection between the st_nearest_points output and nearest reach, but not the point extracted from the st_nearest_points output and the nearest reach.
> pt <- st_set_precision(pt, precision=1000)
> lines <- st_set_precision(lines, precision=1000)
>
> # Select nearest reach
> reach.nearest <- st_nearest_feature(pt, lines)
> reach.nearest <- lines[reach.nearest,]
> st_distance(pt, reach.nearest)
Units: m
[,1]
[1,] 4.684326
>
> # Locate nearest point along nearest reach (Returns line)
> pt.nearest <- st_nearest_points(pt, reach.nearest)
> cat("Line intersects:", st_intersects(pt.nearest, reach.nearest, sparse=FALSE)[1,1], "|")
Line intersects: TRUE |> st_distance(pt.nearest, reach.nearest)
Units: m
[,1]
[1,] 0
>
> pt.nearest <- st_cast(pt.nearest, "POINT")[2]
> cat("Point (st_cast) intersects:", st_intersects(pt.nearest, reach.nearest, sparse=FALSE)[1,1], "|")
Point (st_cast) intersects: FALSE |
> st_distance(pt.nearest, reach.nearest)
Units: m
[,1]
[1,] 0.0001278386I think this problem cannot be solved. Setting precision essentially rounds coordinates to a (scaled) integer grid. If you want to represent points on lines connecting points on an integer grid exactly, you need to store them as rational numbers, i.e. as a ratio of two (possibly very large) integers. What R (or GEOS) does is approximating this ratio when storing it as an eight byte double. By that, we're back at R FAQ 7.31.
bcaradima
commented
6 years ago So if a given point has a nearby location along a line, R cannot exactly match the coordinates of that nearest location even to a limited number of decimal points?
My motive for this particular task is to measure upstream/downstream distances from sites located along a river network, but I suppose I can always work around this by buffering the point slightly and splitting the reach based on that buffer. Interestingly, this can be done in ArcGIS just with intersecting points.
I have since posted an additional question on SE, but I would consider closing it it has no possible solution.
rsbivand
commented
6 years ago Do not use this issues channel, and preferably not SO (sic) either, because your general approach does not appear to have been thought through.
If you are as you now belatedly reveal dealing with points (sampling points) on a stream network, look first at the specific tools designed to that purpose, in particular the SSN package and STARS ArcGIS toolbox (JSS paper at https://doi.org/10.18637/jss.v056.i02), and the openSTARS package using GRASS rather than ArcGIS. None of your interventions so far show that you have actually done enough background work. Do also see Mira Kattwinkel's talk at eRum 2018 on openSTARS.
If having done this and contacted those researchers, you find that precision is not just your misunderstanding, you may be able to contribute to SSN, STARS and/or openSTARS.
bcaradima
commented
6 years ago Thank you for suggesting these interesting alternatives, however there are several reasons why I believe they are unsuitable. I am more interested in quantifying habitat connectivity and quality for fish communities based on their movement ranges, and would like to integrate fish observations, barriers, and other stream attributes as nodes or edges in a directed graph based on a vector-based river network at a large spatial scale. I would prefer to avoid "reinventing the wheel", but I find that these alternatives are either based on closed-source software (STARS) or are not well-suited to the nature of the data (openSTARS) and the types of analyses (SSN) that I'd like to do.
I opened this issue because the st_nearest_points tool returns points that are 10^-11 units from the nearest line and do not intersect the line. @edzer has pointed out that this cannot be resolved if the precision is lowered and the points and lines still do not intersect. My last question (an ignorant one perhaps) before closing would be whether this task can be done at all in R or if it is a limitation of R itself when doing geometric calculations.
Thank you for your helpful replies.
edzer
commented
6 years ago I think it is symptomatic of any system that represents its coordinates with finite-precision floating point numbers rather than ratios of integers not limited by size; see WR Franklin, 1984. Cartographic errors symptomatic of underlying algebra problems; International Symposium on Spatial Data Handling, Zurich, Switzerland 286.
It wouldn't be rocket science to write a package that does this, the question is whether this problem is large enough for you to do so.
> fortunes::fortune("there is no if")
Evelyn Hall: I would like to know how (if) I can extract some of the
information from the summary of my nlme.
Simon Blomberg: This is R. There is no if. Only how.
-- Evelyn Hall and Simon 'Yoda' Blomberg
R-help (April 2005)As @edzer commented, this is a feature of the IEEE 754 floating point representation used on most computing hardware. This is an international standard that R and effectively everyone else chooses to use. Really, your problem is conceptual, as these differences cannot be measured in surveying, and the natural phenomena you are representing with points and lines are not points/lines. If you care to contribute to JTS/GEOS/PostGIS to resolve your problem, fine. Do look at GRASS, which uses different topological choices than JTS/GEOS. There vector points are assigned to streams using a raster stream representation in r.stream.snap, so the rasterization resolution of the stream would decide.
barryrowlingson
commented
6 years ago I suspect the problem is that whether a point lies exactly on a line is defined mathematically, but because of finite precision arithmetic two different algorithms, or even two different implementations of the same algorithm may differ in their opinions. In over-simplified terms, the midpoint of a line is (A+B)/2, but floating point precision might mean that's not the same as (B+A)/2. Then depending on how the code computes the distance from that "midpoint" to the line it might well not be absolute zero.
One algorithm's point on a line is another algorithm's point-not-quite-on-a-line.
bcaradima
commented
6 years ago I've mostly used ArcGIS in the past and locating points along a line was a "trivial" task, with the underlying topology of geodatabases cleanly resolving whether point features intersect. As @edzer and @rsbivand pointed out, differences of 10^-11 units wouldn't matter in my use case, but I just didn't understand why any non-zero distances were returned at all (coming from ArcGIS I haven't seen this happen before).
My current workaround would be to just buffer the sites by a hundredth of a millimeter and then split the lines (which are often hundreds of meters in length) using the buffer. I will also take a look at the possibility of using GRASS.
Again, thank you all for your help.
rsbivand
commented
6 years ago By the way, because ArcGIS is proprietary, you don't see how they address the same issues of computational geometry. They have to make assumptions about snap distances that may or may not be exposed to the user in the programming API.
I have a sites (
pt, POINT) and river network (lines, LINESTRING) dataset. Often, a site is located a few meters away from the nearest line feature. For a given site I would like to locate the nearest point along the nearest line that intersects that line using only thesfpackage (which is a fantastic package BTW!). I have posted this problem on Stack Exchange but have not found a satisfactory and reliable solution since posting. @barryrowlingson and others replied on SE, and @barryrowlingson kindly posted an issue (#788) requesting such a feature onsfthat has already been implemented in 0.6.4.I installed
sf_0.6-4and made use ofst_nearest_feature()andst_nearest_points(), but found that the LINESTRING output from st_nearest_points sometimes does not intersect the nearest line. I'm not sure why this is happening and have tried using thest_snapfunction as well without success. I would appreciate any help in understanding why the nearest point here is still a small distance away from intersecting the nearest line. I've posted a reproducible example below since I'm uncertain if this is a package issue or not; if it isn't then closing the issue would be OK with me.