Scale a lon/lat vector to a certain length in meter

Question

I have some GPS measured lines with 3 points each that are crooked because of uncertainties of the measurement and that I want to replace by lines with a length of exactly 50m centered on the center point of the measured lines and with the general direction of the measured points.

To do this, I want to go 25m from the center point (given by lon/ lat coordinates) into a certain direction, which is given by a lon/ lat vector. Now the problem is, that at different latitudes 25m are not the same in degrees. I did some calculation to get the length in degrees that 1m is at the location:

# 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

I have a dataframe (coords) which contains in each row the coordinates of the start s, center c and end e point of a line. I did a regression to get the line that approximates the points' general direction best (here, I used [1] for the first row of the data; if it works it will obviously be looped to repeat the process for all the lines):

# 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]

Now I tried to go in the direction given by the model and scale the vector by whatever it takes to make it 25m long at the given latitude, but somehow the new start and end points are not on the regression line anymore:

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

Does someone see where the error comes in or know how to resolve the problem? I guess I didn't scale the vector properly...

In case you want the first line of the data to test this:

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

Not sure i what format you want the output, but this might get you started...

results in two new columns; new_s and new_e, in lon,lat - format..

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 ) 
                        )
  )

output on a map

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" )

old points are red, new points are blue.. center point stays the same....