2

I have been trying to git this right, in fact it is working but the distance is wrong:

  • 1- I have a collection of points
  • 2- I have a point within which I need to find all points that are 200m far
  • 3- Normally, It should return 6 points, however it returns only 4
library(sf)
library(geosphere)

​
#My points
x<-c(2.3009132,2.2999853,2.2995872,2.2985374,2.2991502,2.2984043,2.3054471,2.3009132,2.3048155,2.3014964)
y<-c(48.8511847,48.8505062,48.8502346,48.8495305,48.8499405,48.8494376,48.8542721,48.8511847,48.853842,48.8515819)​
y<-c(48.8511847,48.8505062,48.8502346,48.8495305,48.8499405,48.8494376,48.8542721,48.8511847,48.853842,48.8515819)
df<-data.frame(x=x,y=y)
#Transforming to SF object 
sdf<-st_transform(st_as_sf(df, coords = c("x", "y"), 
                           crs = 4326, agr = "constant"),3857)
#My point to which I need to calculte
pnt<- st_transform( 
    st_sfc(st_point(x = c(2.3009132, 48.8511847)), crs = 4326), 3857)
#A buffer of 200m arround my point
buffer <- st_buffer(pnt,200)
#getting points within the buffer
intr <- st_intersects(sdf, buffer, sparse=F)
#transforming back to lon/lat
sdf <- st_transform(sdf,4326)

#getting the selected points
sdf<-sdf[which(unlist(intr)),]

#Only 4 points were found
> sdf
Simple feature collection with 4 features and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 256033.2 ymin: 6249533 xmax: 256201.4 ymax: 6249715
epsg (SRID):    3857
proj4string:    +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
                  geometry
1 POINT (256136.5 6249648)
2 POINT (256033.2 6249533)
3 POINT (256136.5 6249648)
4 POINT (256201.4 6249715)

#To verify I have calculated the distance to my point
t_sdf<-df%>% mutate(d = by(df, 1:nrow(df), function(row) {
+     distHaversine(c(row$x, row$y), c(2.3009132, 48.8511847), r = 6378137)
+ }))
#6 points are less than 200m to my point
> t_sdf %>% arrange(d)
          x        y         d
1  2.300913 48.85118   0.00000
2  2.300913 48.85118   0.00000
3  2.301496 48.85158  61.48172
4  2.299985 48.85051 101.61024
5  2.299587 48.85023 143.59844
6  2.299150 48.84994 189.36954
.....
roger
  • 89
  • 10

2 Answers2

3

Straight line (Euclidean) distances between a certain set of points vary across different projections.

As discussed here, 3857 is not a great choice for distance calculations. If we use a projection built for Northern France, where the points are located (perhaps EPSG:27561 - NTF (Paris) / Lambert Nord France), we get the expected results.

sdf <- st_as_sf(df, coords = c("x", "y"), 
                crs = 4326, agr = "constant")
sdf_proj <- st_transform(sdf, 3857)
sdf_France_proj <- st_transform(sdf, 27561)

Great circle distances

> st_distance(pnt, sdf, by_element = T) %>% sort()
Units: [m]
 [1]   0.00000   0.00000  61.50611 101.64007 143.64481 189.43485 253.46142 267.67954 411.50933 478.11306

Mercator (3857) Euclidean distances

> st_distance(pnt_proj, sdf_proj, by_element = T) %>% sort()
Units: [m]
 [1]   0.00000   0.00000  93.43522 154.41781 218.22699 287.78463 385.04306 406.64205 625.14635 726.33083

France projection (27561) Euclidean distances

> st_distance(pnt_France_proj, sdf_France_proj, by_element = T) %>% sort()
Units: [m]
 [1]   0.00000   0.00000  61.50289 101.63477 143.63731 189.42497 253.44822 267.66559 411.48772 478.08794

And here are the plots with buffers:

library(gridExtra)

buffer_proj <- st_buffer(pnt_proj, 200)
buffer_France_proj <- st_buffer(pnt_France_proj, 200)

proj <- ggplot() + g
eom_sf(data = buffer_proj, fill = "blue", alpha = 0.5) + 
  geom_sf(data = sdf_proj, color = "red") + 
  ggtitle("Mercator (3857)")

France_proj <- ggplot() + 
  geom_sf(data = buffer_France_proj, fill = "blue", alpha = 0.5) + 
  geom_sf(data = sdf_proj, color = "red") + 
  ggtitle("France proj (27561)")

grid.arrange(proj, France_proj, ncol = 1)

enter image description here

Eugene Chong
  • 1,591
  • 5
  • 16
  • Thank you, this was exactly what it caused this problem – roger Dec 13 '19 at 16:32
  • Please, for the points I got, If I want to add a column that is the distance of each point to the first point in the box(because its a street, it should be the first point in the street) – roger Dec 13 '19 at 16:44
  • 1
    We actually discussed this in a previous post of yours. See the bottom of my answer here - look at mutate(): https://stackoverflow.com/questions/59250592/distance-between-all-the-points-to-the-first-point/59250811#59250811 – Eugene Chong Dec 13 '19 at 16:46
  • Yes you are right, however, there are thousands of points, I'm ordering everytime the points then I calculate the distance to the first. I'm new in geomatics, I know I can access to xmin, ymin by using st_bbox Is the first point(in a given Street) always c(xmin,ymin) (So I will no longer order to calculate)? – roger Dec 13 '19 at 16:54
  • I think you should not need to order the points first. All you need to do is find the first one on every street, find its distance to all the other points, and then sort by distance. But I may be misunderstanding your problem. If that's true, I recommend you create a new post and explain in further detail. – Eugene Chong Dec 13 '19 at 17:02
  • In fact, I have multiple streets, is there an automatical way to get the first point? In fact in all the sf objects, the bounding boxes are accessible via: st_bbox(df) – roger Dec 13 '19 at 17:08
  • 1
    that is possibly a solution, but without seeing your data more closely, it's hard to say. I think you should create a new post that includes more of your actual data. This seems too complicated to discuss in the comments, and it'll get a bigger audience. – Eugene Chong Dec 13 '19 at 17:11
2

I don't know how you are calculating the distance to your point. But when when you do the calculations with the data you provided there are only 4 points that fall into the circle. Maybe you can use QGIS or ArcGIS to verify this. But here is how your 4 points look like within the circle:

enter image description here

Aliyan Haq
  • 302
  • 1
  • 7
  • distHaversine Calculates he shortest distance between two points (i.e., the 'great-circle-distance' or 'as the crow flies'), I have also tested distGeo and other they all give the same result. – roger Dec 13 '19 at 15:08
  • 2
    FWIW, I've taken this data and calculated straight line distance as well using a separate implementation of haversine distance and get the same result as OP in the question. Incidentally, I also get the same result as this answer when viewing the data. – ander2ed Dec 13 '19 at 15:12
  • 1
    Very weird. One thing to note though is that when I plotted the data. Even though everything was in meters (units=m), the x and y axis were in lat long. I have edited my picture above. – Aliyan Haq Dec 13 '19 at 15:26
  • 1
    Maybe try using a different projection – Aliyan Haq Dec 13 '19 at 16:15