将 lon/lat 矢量缩放到一定长度(以米为单位)
Scale a lon/lat vector to a certain length in meter
我有一些 GPS 测量线,每条线有 3 个点,由于测量的不确定性,我想用长度正好为 50m 的线来替换测量线的中心点和测量点的大致方向。
为此,我想从中心点(由 lon/lat 坐标给出)向 25m 移动到某个方向,该方向由 lon/lat 向量给出.现在的问题是,在不同的纬度 25m 度数是不一样的。我做了一些计算以获得 1m 在该位置的度数长度:
# 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
我有一个 dataframe (coords),它在每一行中包含开始的坐标 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]
现在我尝试按照模型给定的方向前进,并根据给定纬度使矢量 25m 长所需的任何比例缩放矢量,但不知何故新 start和end点不再在回归线上:
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")
不确定你想要什么格式的输出,但这可能会让你开始......
产生两个新列; new_s
和 new_e
,采用经度、纬度格式..
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" )
旧点为红色,新点为蓝色..中心点保持不变....
我有一些 GPS 测量线,每条线有 3 个点,由于测量的不确定性,我想用长度正好为 50m 的线来替换测量线的中心点和测量点的大致方向。
为此,我想从中心点(由 lon/lat 坐标给出)向 25m 移动到某个方向,该方向由 lon/lat 向量给出.现在的问题是,在不同的纬度 25m 度数是不一样的。我做了一些计算以获得 1m 在该位置的度数长度:
# 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
我有一个 dataframe (coords),它在每一行中包含开始的坐标 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]
现在我尝试按照模型给定的方向前进,并根据给定纬度使矢量 25m 长所需的任何比例缩放矢量,但不知何故新 start和end点不再在回归线上:
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")
不确定你想要什么格式的输出,但这可能会让你开始......
产生两个新列; new_s
和 new_e
,采用经度、纬度格式..
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" )
旧点为红色,新点为蓝色..中心点保持不变....