How to find the row number based on a set of SFC coordinates, in order to pull out the value from another column

Question

I am new to geodetic data in R, so apologies if this is way off. I have got a dataset, with a column geometry, which contains SFC points, and a second column IntersectingPoints:

IntersectingPoints             geometry

<list>                      <POINT [°]>
<int [4,415]>      (144.9179 -45.59569)
<int [3,388]>     (5.740768 49.15625)
<int [6,214]>        (2.18701 52.80779)
<int [7,948]>       (3.108628 50.58135)
<int [2,298]>      (-96.77081 37.37828)

I would like to input an x and y value, in order to find the row number, which I can then use to pull out the same value from the IntersectingPoints column. The code I have written is:

x = 144.9179
y = -45.59569
point_ref=st_point(c(x,y))
which(pt$geometry == point_ref)

I was expecting this to output the first row, however, this does not work, I think as the geometry is classed as a list.

How I can output the row number?

EDIT

The error when trying to run my code was:

Error in which(pt$geometry == point_ref) : 
  'list' object cannot be coerced to type 'double'
In addition: Warning messages:
1: In which(pt$geometry == point_ref) :
  Incompatible methods ("Ops.sfc", "Ops.sfg") for "=="
2: In pt$geometry == point_ref :
  longer object length is not a multiple of shorter object length

When I run st_sfc(point_ref) == st_geometry(pt), this does not find any results.

Output for 'st_sfc(point_ref):

Geometry set for 1 feature
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 144.9179 ymin: -45.59569xmax: 144.9179 ymax: -45.59569
CRS:           NA
POINT (144.9179 -45.59569)

Output for 'st_geometry(pt)' is:

Geometry set for 139709 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -172.9813 ymin: -56.81466 xmax: 171.9221 ymax: 71.08403
Geodetic CRS:  WGS 84
First 5 geometries:
POINT (144.9179 -45.59569)

I can see that the points are the same, so I would have thought that when these were run together, a match would have been found!

Answer 1

You are comparing an st_point object (which is an sfg object) with a vector of st_point objects wrapped in the st_sfc class, which holds metadata about coordinate systems etc. You don't say, but I guess you are getting this error:

> point_ref == st_geometry(pt)
Error: 'list' object cannot be coerced to type 'double'
In addition: Warning message:
Incompatible methods ("Ops.sfg", "Ops.sfc") for "==" 

You can turn your point_ref into an sfc and then compare:

> st_sfc(point_ref) == st_geometry(pt)
 [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

and that gives a standard R boolean you can use to extract from rows of pt, or run which on it to get row numbers etc.

HOWEVER....

Doing this is comparing floating point numbers, which is often an ill-advised thing to do. Depending on where your numbers came from, numbers that look equal, or even mathematically should be equal aren't. eg:

> sqrt(2) * sqrt(2) == 2
[1] FALSE

An example using points would be this. x and xx look the same:

> x
[1] 144.9179
> xx
[1] 144.9179

Comparing a point made with x and y with another point also made with x and y gives equality:

> st_point(c(x,y)) == st_point(c(x,y))
[1] TRUE

but with x and xx fails equality:

> st_point(c(x,y)) == st_point(c(xx,y))
[1] FALSE

even though these two things look the same:

> st_point(c(x,y))
POINT (144.9179 -45.59569)
> st_point(c(xx,y))
POINT (144.9179 -45.59569)

The difference is tiny:

> xx - x
[1] 9.998757e-11

but it will cause equality tests to fail.

Unless your test point is an unmodified subset of your comparison set, you can't rely on equality testing. Example:

Construct 26 points with random coords:

> d = st_as_sf(data.frame(x=runif(26), y=runif(26), ID=LETTERS), coords=1:2)

Take the 12th point and compare it to all the points. Twelfth value returned is TRUE:

> p = st_geometry(d)[12]
> p == st_geometry(d)
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE FALSE

But if I look at that point and make a point with those numbers:

> p
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 0.2990404 ymin: 0.005280197 xmax: 0.2990404 ymax: 0.005280197
CRS:           NA
POINT (0.2990404 0.005280197)
> pfail = st_sfc(st_point(c(0.2990404,0.005280197)))

it doesn't match:

> pfail == st_geometry(d)
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE FALSE

The distance is:

> min(st_distance(pfail, st_geometry(d))[1,])
[1] 1.727231e-08

so if you wanted to test if two points were "the same" you need to test within some threshold to be safe against these kinds of approximations.

See R FAQ 7.31 for more details.

https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f