Определить точки на указанном расстоянии в R

Question

Из списка из 10 000 станций с десятичными координатами я пытаюсь идентифицировать станции, находящиеся на расстоянии 100 футов друг от друга, исходя из расстояния, вычисленного между этими станциями, и создания подмножества этих станций. В последнем списке я хочу иметь имена станций, которые находятся на расстоянии 100 футов друг от друга, их широты и долготы и расстояния между ними.

Я нашел аналогичные вопросы для других платформ, таких как MathWorks (с использованием rangesearch) или в SQL или JAVA, но ни один из R.

Есть ли способ сделать это в R? Самый близкий ответ, который я нашел, был в Listing number of obervations by location, который перечисляет количество наблюдений на расстоянии, но кажется, что ответы были неполными и не могут определить станции, находящиеся на определенном расстоянии друг от друга.

В основном я пытаюсь выяснить, какие станции расположены.

Я бы очень признателен за любую помощь в этом.

Answer 1

два подхода.

Первый создает матрицу расстояний, используя earth.dist(...) в пакете fossil, а затем использует data.tables для сборки таблицы результатов.

Второй использует distHaversine(...) в пакете geosphere для расчета расстояний и сборки окончательной таблицы колокаций за один шаг. Последний подход может быть или не быть быстрее, но, безусловно, будет более эффективным с точки зрения памяти, поскольку он никогда не хранит полную матрицу расстояний. Кроме того, этот подход поддается использованию других дистанционных мер в geosphere, например, distVincentySphere(...), distVincentyEllipsoid(...) или distMeeus(...).

Обратите внимание, что фактические расстояния немного отличаются, вероятно, потому, что earth.dist(...) и distHaversine(...) используют несколько разные оценки радиуса земли. Также обратите внимание, что оба подхода здесь основаны на номерах станций для идентификаторов. Если станции имеют имена, код должен быть слегка изменен.

Первый подход: Использование earth.dist(...)

df = read.table(header=T,text="long lat 
       1 -74.20139 39.82806 
       2 -74.20194 39.82806 
       3 -74.20167 39.82806 
       4 -74.20197 39.82824 
       5 -74.20150 39.82814 
       6 -74.26472 39.66639 
       7 -74.17389 39.87111 
       8 -74.07224 39.97353 
       9 -74.07978 39.94554")    # your sample data 
library(fossil)          # for earth.dist(...) 
library(data.table) 
sep.ft <- 200          # critical separation (feet) 
sep.km <- sep.ft*0.0003048      # critical separation (km) 
m  <- as.matrix(earth.dist(df))    # distance matrix in km 
coloc <- data.table(which(m<sep.km, arr.ind=T)) # pairs of stations with dist<200 ft 
setnames(coloc,c("row","col"),c("ST.1","ST.2"))  # rename columns to reflect station IDs 
coloc <- coloc[ST.1<ST.2,]      # want only lower triagular part 
coloc[,dist:=m[ST.1,ST.2]/0.0003048,by="ST.1,ST.2"] # append distances in feet 
remove(m)           # don't need distance matrix anymore... 
stations <- data.table(id=as.integer(rownames(df)),df) 
setkey(stations,id) 
setkey(coloc,ST.1) 
coloc[stations,c("long.1","lat.1"):=list(long,lat),nomatch=0] 
setkey(coloc,ST.2) 
coloc[stations,c("long.2","lat.2"):=list(long,lat),nomatch=0] 

Производит это:

coloc 
#  ST.1 ST.2  dist long.1 lat.1 long.2 lat.2 
# 1: 1 2 154.13436 -74.20139 39.82806 -74.20194 39.82806 
# 2: 1 3 78.46840 -74.20139 39.82806 -74.20167 39.82806 
# 3: 2 3 75.66596 -74.20194 39.82806 -74.20167 39.82806 
# 4: 1 4 175.31180 -74.20139 39.82806 -74.20197 39.82824 
# 5: 2 4 66.22069 -74.20194 39.82806 -74.20197 39.82824 
# 6: 3 4 106.69018 -74.20167 39.82806 -74.20197 39.82824 
# 7: 1 5 42.45634 -74.20139 39.82806 -74.20150 39.82814 
# 8: 2 5 126.71608 -74.20194 39.82806 -74.20150 39.82814 
# 9: 3 5 55.87449 -74.20167 39.82806 -74.20150 39.82814 
# 10: 4 5 136.67612 -74.20197 39.82824 -74.20150 39.82814 

Второй подход: Использование distHaversine(...)

library(data.table) 
library(geosphere) 
sep.ft <- 200      # critical separation (feet) 
stations <- data.table(id=as.integer(rownames(df)),df) 

d <- function(x){      # distance between station[i] and all subsequent stations 
    r.ft <- 6378137*3.28084    # radius of the earth, in feet 
    if (x[1]==nrow(stations)) return() # don't process last row 
    ref <- stations[(x[1]+1):nrow(stations),] 
    z <- distHaversine(ref[,2:3,with=F],x[2:3], r=r.ft) 
    z <- data.table(ST.1=x[1], ST.2=ref$id, dist=z, long.1=x[2], lat.1=x[3], long.2=ref$long, lat.2=ref$lat) 
    return(z[z$dist<sep.ft,]) 
} 
coloc.2 = do.call(rbind,apply(stations,1,d)) 

Производит это:

coloc.2 
#  ST.1 ST.2  dist long.1 lat.1 long.2 lat.2 
# 1: 1 2 154.26350 -74.20139 39.82806 -74.20194 39.82806 
# 2: 1 3 78.53414 -74.20139 39.82806 -74.20167 39.82806 
# 3: 1 4 175.45868 -74.20139 39.82806 -74.20197 39.82824 
# 4: 1 5 42.49191 -74.20139 39.82806 -74.20150 39.82814 
# 5: 2 3 75.72935 -74.20194 39.82806 -74.20167 39.82806 
# 6: 2 4 66.27617 -74.20194 39.82806 -74.20197 39.82824 
# 7: 2 5 126.82225 -74.20194 39.82806 -74.20150 39.82814 
# 8: 3 4 106.77957 -74.20167 39.82806 -74.20197 39.82824 
# 9: 3 5 55.92131 -74.20167 39.82806 -74.20150 39.82814 
# 10: 4 5 136.79063 -74.20197 39.82824 -74.20150 39.82814 

Answer 2

Вот лишь некоторые случайные примеры данных

set.seed(1234) 
x= sample(1:100,50) 
y= sample(1:100,50) 
M=cbind(x,y) 
plot(M) 

Вы вычислить расстояния в виде матрицы таким образом, что исходные строки могут быть извлечены легко. Это может быть сделано с помощью функции which с arr.ind = T, например, так:

DM= as.matrix(dist(M)) 
neighbors=which(DM < 5, arr.ind=T) 
neighbors= neighbors[neighbors[,1]!=neighbors[,2]] 

Таким образом, вы можете определить точки, которые говорят менее 5 единиц евклидов расстояния (отдельно после удаления собственных самостоятельных отношений):

points(M[neighbors,], col="red")