Масштабируйте вектор долгота / широта до определенной длины в метрах

Question

У меня есть несколько линий, измеренных GPS, с 3 точками каждая, которые искривлены из-за неточностей измерения и которые я хочу заменить линиями длиной точно 50 м с центром в центральной точке измеряемого объекта. линий и с общим направлением измеренных точек.

Для этого я хочу go 25 м от центральной точки (заданной координатами долгота / широта) в определенном направлении , который задается вектором долгота / широта. Теперь проблема в том, что на разных широтах 25 м не совпадают в градусах. Я сделал некоторые вычисления, чтобы получить длину в градусах, в которой 1 м находится в месте:

# Calculate the length of a 1 degree difference
earth_circumference <- 40075016.7
for(i in 1:nrow(coords)){
  coords$d_m_lon[i] <- earth_circumference/360*cos(abs(coords$c_lat[1])*pi/180)
  coords$d_m_lat[i] <- earth_circumference/360
}

# Claculate the difference in degree that 1m at the surface makes
coords$m_d_lon <- 1/coords$d_m_lon
coords$m_d_lat <- 1/coords$d_m_lat

У меня есть фрейм данных ( координаты ), который содержится в каждой строке координаты начала s , центра c и конца e точки линии. Я сделал регрессию, чтобы получить линию, которая наилучшим образом аппроксимирует общее направление точек (здесь я использовал [1] для первой строки данных; если это сработает, очевидно, что цикл будет повторяться, чтобы повторить процесс для всех строк):

# Calculate average transect directions
points <- data.frame(
  x = c(coords$s_lon[1], coords$c_lon[1], coords$e_lon[1]),
  y = c(coords$s_lat[1], coords$c_lat[1], coords$e_lat[1])
)
mod <- lm(points$y ~ points$x)
a <- coefficients(mod)[1]
b <- coefficients(mod)[2]

Теперь я попытался изменить go в направлении, заданном моделью, и масштабировать вектор, чтобы сделать его 25 м в длину заданная широта, но каким-то образом новые точки начало и конец больше не находятся на линии регрессии:

center <- c(coords$c_lon[1], coords$c_lat[1])
direction <- c(b-a, 1)/sqrt((b-a)^2 + 1)

start <- center + 25*c(coords$m_d_lon[1], coords$m_d_lat[1])*direction
end <- center - 25*c(coords$m_d_lon[1], coords$m_d_lat[1])*direction

Кто-нибудь видит, где возникает ошибка, или знает как решить проблему? Думаю, я неправильно масштабировал вектор ...

Если вы хотите, чтобы первая строка данных проверяла это:

coords <- structure(list(s_lat = -29.6032, s_lon = 29.3376, c_lat = -29.6032, 
    c_lon = 29.3379, e_lat = -29.6032, e_lon = 29.3381, d_m_lon = 96788.6617220582, 
    d_m_lat = 111319.490833333, m_d_lon = 1.03317886848321e-05, 
    m_d_lat = 8.98315283796251e-06), row.names = 1L, class = "data.frame")

Answer 1

Не уверен, какой формат вы хотите получить, но это может помочь вам начать ...

приводит к появлению двух новых столбцов; new_s и new_e, в формате lon, lat ..

library( tidyverse )  
library( geosphere )

coords <- coords %>%
  #calculate bearing from center ppint to s-point
  mutate( bearing_c_to_s = pmap( list ( a = c_lon,
                                        b = c_lat,
                                        x = s_lon,
                                        y = s_lat ),
                                 ~ geosphere::bearing( c(..1, ..2), c(..3, ..4) ) 
                                 )
          ) %>%
  #calculate bearing from center ppint to s-point
  mutate( bearing_c_to_e = pmap( list ( a = c_lon,
                                        b = c_lat,
                                        x = e_lon,
                                        y = e_lat ),
                                 ~ geosphere::bearing( c(..1, ..2), c(..3, ..4) )  
                                 )
          ) %>%
  #calculate new point s coordinates, 25 meters from c using bearing c_to_s
  mutate( new_s = pmap( list ( a = c_lon,
                               b = c_lat,
                               x = bearing_c_to_s,
                               y = 25 ),
                        ~ geosphere::destPoint( c(..1, ..2), ..3, ..4 ) 
                        )
          ) %>%
  #calculate new point e coordinates, 25 meters from c using bearing c_to_e
  mutate( new_e = pmap( list ( a = c_lon,
                               b = c_lat,
                               x = bearing_c_to_e,
                               y = 25 ),
                        ~ geosphere::destPoint( c(..1, ..2), ..3, ..4 ) 
                        )
  )

вывод на карту

library( sf )
library( data.table )
setDT(coords)
coords_old <- data.table::melt(coords, measure.vars = patterns( lat = "^[a-z]_lat", lon = "^[a-z]_lon" ) )[, c("lon","lat")]
coords_old <- sf::st_as_sf( coords_old, coords = c("lon", "lat"), crs = 4326 )
coords[, c("new_s_lon", "new_s_lat") := ( as.list( unlist( new_s ) ) ) ]
coords[, c("new_e_lon", "new_e_lat") := ( as.list( unlist( new_e ) ) ) ]
coords_new <- data.table::melt(coords, measure.vars = patterns( lat = "^(c|new_[se])_lat", lon = "^(c|new_[se])_lon" ) )[, c("lon","lat")]
coords_new <- sf::st_as_sf( coords_new, coords = c("lon", "lat"), crs = 4326 )

library( leaflet )
leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers( data = coords_old, color = "red" ) %>%
  addCircleMarkers( data = coords_new, color = "blue" )

старые точки красные, новые точки синие .. центральная точка остается прежней ....