距点一定距离内的点
Points within a distance from a point
我一直在努力git这个对,其实是可以的,但是距离不对:
- 1- 我有 collection 分
- 2- 我有一个点,我需要在其中找到所有点
200 米远
- 3- 通常情况下,它应该 return 6 点,但它 return 只有 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
.....
我不知道你是如何计算到你的点的距离的。但是当你用你提供的数据进行计算时,只有 4 个点落入圆圈中。也许你可以使用 QGIS 或 ArcGIS 来验证这一点。但这是你的 4 个点在圆圈内的样子:
一组特定点之间的直线(欧几里德)距离因不同的投影而异。
正如所讨论的 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),我们得到了预期的结果。
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)
大圆距离
> 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
墨卡托 (3857) 欧氏距离
> 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
法国投影 (27561) 欧氏距离
> 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
这里是带缓冲区的图:
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)
我一直在努力git这个对,其实是可以的,但是距离不对:
- 1- 我有 collection 分
- 2- 我有一个点,我需要在其中找到所有点 200 米远
- 3- 通常情况下,它应该 return 6 点,但它 return 只有 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
.....
我不知道你是如何计算到你的点的距离的。但是当你用你提供的数据进行计算时,只有 4 个点落入圆圈中。也许你可以使用 QGIS 或 ArcGIS 来验证这一点。但这是你的 4 个点在圆圈内的样子:
一组特定点之间的直线(欧几里德)距离因不同的投影而异。
正如所讨论的 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),我们得到了预期的结果。
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)
大圆距离
> 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
墨卡托 (3857) 欧氏距离
> 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
法国投影 (27561) 欧氏距离
> 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
这里是带缓冲区的图:
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)