Найти ближайшее расстояние от пространственной точки с указанным направлением

Question

Я хотел бы рассчитать ближайшее расстояние от пространственной точки до пространственных линий (или многоугольников) для заданных ориентиров (0,45,90,135,180,225,270,315).

Идея состоит в том, чтобы рассчитать индекс воздействия для числазаливов вдоль береговой линии.Ниже приведен простой пример:

Создание линий

library(sp)
coords<-structure(list(lon = c(-6.1468506, -3.7628174, -3.24646, 
-3.9605713, -4.4549561, -4.7955322, -4.553833, -5.9710693, -6.1468506), 
lat = c(53.884916, 54.807017, 53.46189, 53.363665, 53.507651, 53.363665, 53.126998, 53.298056,53.884916)), class = "data.frame", row.names = c(NA,-9L))
l<-Line(coords)
sl<-SpatialLines(list(Lines(list(l),ID="a")),proj4string=CRS("+init=epsg:4326"))

Создание точки

pt<-SpatialPoints(coords[5,]+0.02,proj4string=CRS("+init=epsg:4326"))

Участок

plot(sl)
plot(pt,add=T)

У меня проблемынайти примеры того, каким может быть следующий шаг, и нужна помощь. Пример того, какое расстояние я хотел бы рассчитать

Answer 1

В качестве справки: я не герой, когда дело доходит до геоданных в R.

Также: я не автоматизировал расчет для всех подшипников, но выполнял операции вручную, чтобы получить расстояниечтобы пересечь на 45-й подшипник.
Вам придется выяснить петли самостоятельно, так как у меня нет времени.Не стесняйтесь предоставить / опубликовать свои окончательные результаты / код здесь, когда вы закончите.

Вот мой шаг в решении этой проблемы, шаг за шагом.

#load libraries used
library(geosphere)
library(tidyverse)
library(sf)

#get bearings of lines of the polygon
df.poly <- coords %>%
  mutate( lon_next = lead(lon), lat_next = lead(lat) ) %>%
  mutate( bearing_to_next = ifelse( !is.na( lon_next ), 
                                    unlist( pmap( list( a = lon, b = lat, x = lon_next, y = lat_next ),
                                                  ~ round( geosphere::bearing( c(..1, ..2), c(..3, ..4) ) )
                                                  )
                                            ),
                                    NA ) 
          ) %>%
  filter( !is.na( lon_next ) )

#         lon      lat bearing_to_next
# 1 -6.146851 53.88492              56
# 2 -3.762817 54.80702             167
# 3 -3.246460 53.46189            -103
# 4 -3.960571 53.36366             -64
# 5 -4.454956 53.50765            -125
# 6 -4.795532 53.36366             148
# 7 -4.553833 53.12700             -78
# 8 -5.971069 53.29806             -10

#find intersection point based on the intersection of two 'great circles' 
#from two points with a bearing
gcIntersectBearing( 
  #coordinates 2nd point of polyline, with bearing to third point
  c( -3.7628174, 54.807017 ), 167,  
  #coordinates point, with bearing of 45
  c( -4.454956, 53.50765 ), 45 )

#            lon      lat      lon       lat
# [1,] -3.476074 54.07798 176.5239 -54.07798

давайте посмотрим, что мыдошли до сих пор

p_intersect <- data.frame( lon = -3.476074, lat = 54.07798 ) %>% 
  st_as_sf( coords = c( "lon", "lat" ), crs = 4326 )

startpoint <- coords %>% slice(5) %>% mutate( lon = lon + 0.02, lat = lat + 0.02 ) %>%
  st_as_sf( coords = c("lon","lat"), crs = 4326 )

poly <- coords %>%
  as.matrix() %>%
  list() %>%
  st_polygon() %>%
  st_sfc() %>%
  st_set_crs( 4326 )

mapview::mapview( list(poly, startpoint, p_intersect) )

Расположение точки пересечения p_intersect на многоугольнике poly от startpoint с45-градусный подшипник выглядит правильно.

Теперь вы можете рассчитать расстояние следующим образом:

#calculate distance
st_distance( startpoint, p_intersect )
# Units: [m]
#         [,1]
# [1,] 87993.3

Карты Google, похоже, согласуются с расстоянием (небольшая разница из-за щелчка мышью по точкам, но для меня это нормально)

Теперь вам нужно придумать какую-нибудь хитрую петлю / векторизацию, и все готово :) Я должен вернуться к своей настоящей работе.

Answer 2

Вы можете использовать библиотеку geosphere для этого. Вам нужно будет добавить CRS к вашим очкам, хотя:

library(geosphere)

pt <- SpatialPoints(c[5,],
                    proj4string=CRS("+init=epsg:4326")) 

А затем используйте dist2Line функцию:

st_distance(st_cast(sl, "POINT"), pt)

#     distance       lon      lat ID
#[1,] 2580.843 -4.451901 53.50677  1

---

В качестве альтернативы вы можете преобразовать свои полилинии в точки, используя пакет sf, а затем получить матрицу расстояний (вам нужно преобразовать ваши объекты в класс sf):

library(sf)

sl <- SpatialLines(list(Lines(list(l),ID="a")),
                   proj4string=CRS("+init=epsg:4326")) %>% 
  st_as_sf()

pt <- SpatialPoints(coords[5,]+0.02,
                    proj4string=CRS("+init=epsg:4326")) %>% 
  st_as_sf()

st_distance(st_cast(sl, "POINT"), pt)

#Units: [m]
#            [,1]
# [1,] 119833.165
# [2,] 149014.814
# [3,]  79215.071
# [4,]  36422.390
# [5,]   2591.267
# [6,]  30117.701
# [7,]  45287.637
# [8,] 105289.230
# [9,] 119833.165

Answer 3

Спасибо @patL и @Wimpel, я использовал ваши предложения, чтобы найти решение этой проблемы.

Сначала я создаю пространственные линии заданного расстояния и пеленгов от начальной точки, используя destPoint::geosphere. Затем я использую gIntersection::rgeos, чтобы получить пространственные точки, где каждый разрез пересекает береговую линию. Наконец, я вычисляю расстояние от исходной точки до всех точек пересечения для каждой линии пересечения, соответственно, используя gDistance::rgeos, и устанавливаю минимальное значение, т.е. ближайшее пересечение.

загрузка пакетов

pkgs=c("sp","rgeos","geosphere","rgdal") # list packages
lapply(pkgs,require,character.only=T) # load packages

создать данные

побережье

coords<-structure(list(lon =c(-6.1468506,-3.7628174,-3.24646, 
-3.9605713,-4.4549561,-4.7955322,-4.553833,-5.9710693,-6.1468506), 
lat=c(53.884916,54.807017,53.46189,53.363665,53.507651,53.363665,53.126998,53.298056,53.884916)), class = "data.frame", row.names = c(NA,-9L))
l=Line(coords)
sl=SpatialLines(list(Lines(list(l),ID="a")),proj4string=CRS("+init=epsg:4326"))

точка

sp=SpatialPoints(coords[5,]+0.02,proj4string=CRS("+init=epsg:4326"))
p=coordinates(sp) # needed for destPoint::geosphere

создать линии разреза

b=seq(0,315,45) # list bearings

tr=list() # container for transect lines

for(i in 1:length(b)){
    tr[[i]]<-SpatialLines(list(Lines(list(Line(list(rbind(p,destPoint(p=p,b=b[i],d=200000))))),ID="a")),proj4string=CRS("+init=epsg:4326")) # create spatial lines 200km to bearing i from origin
    }

расчет расстояния

minDistance=list() # container for distances


for(j in 1:length(tr)){ # for transect i
    intersects=gIntersection(sl,tr[[j]]) # intersect with coastline
    minDistance[[j]]=min(distGeo(sp,intersects)) # calculate distances and use minimum
}

do.call(rbind,minDistance)

В действительности исходной точкой является фрейм данных пространственной точки, и этот процесс повторяется многократно для нескольких сайтов. При выполнении пересечения также имеется ряд результатов NULL, поэтому цикл включает в себя оператор if.