График круга с определенным радиусом вокруг точки на карте в ggplot2
У меня есть карта с 8 пунктами, нанесенными на нее:
library(ggplot2)
library(ggmap)
data = data.frame(
ID = as.numeric(c(1:8)),
longitude = as.numeric(c(-63.27462, -63.26499, -63.25658, -63.2519, -63.2311, -63.2175, -63.23623, -63.25958)),
latitude = as.numeric(c(17.6328, 17.64614, 17.64755, 17.64632, 17.64888, 17.63113, 17.61252, 17.62463))
)
island = get_map(location = c(lon = -63.247593, lat = 17.631598), zoom = 13, maptype = "satellite")
islandMap = ggmap(island, extent = "panel", legend = "bottomright")
RL = geom_point(aes(x = longitude, y = latitude), data = data, color = "#ff0000")
islandMap + RL + scale_x_continuous(limits = c(-63.280, -63.21), expand = c(0, 0)) + scale_y_continuous(limits = c(17.605, 17.66), expand = c(0, 0))
Теперь я хочу построить круг вокруг каждого из 8 построенных по графику местоположений. Круг должен иметь радиус 450 метров.
Это то, что я имею в виду, но затем используя ggplot: https://gis.stackexchange.com/info/119736/ggmap-create-circle-symbol-where-radius-represents-distance-miles-or-km
Как я могу это достичь?
Ответы
Ответ 1
Если вы работаете только на небольшой площади земли, вот приблизительное. Каждый градус широты составляет 40075/360 километров. Каждая градус долготы представляет (40075/360) * cos (широта) километр. При этом мы можем рассчитать приблизительно кадр данных, включающий все точки на кругах, зная центры окружности и радиус.
library(ggplot2)
library(ggmap)
data = data.frame(
ID = as.numeric(c(1:8)),
longitude = as.numeric(c(-63.27462, -63.26499, -63.25658, -63.2519, -63.2311, -63.2175, -63.23623, -63.25958)),
latitude = as.numeric(c(17.6328, 17.64614, 17.64755, 17.64632, 17.64888, 17.63113, 17.61252, 17.62463))
)
#################################################################################
# create circles data frame from the centers data frame
make_circles <- function(centers, radius, nPoints = 100){
# centers: the data frame of centers with ID
# radius: radius measured in kilometer
#
meanLat <- mean(centers$latitude)
# length per longitude changes with lattitude, so need correction
radiusLon <- radius /111 / cos(meanLat/57.3)
radiusLat <- radius / 111
circleDF <- data.frame(ID = rep(centers$ID, each = nPoints))
angle <- seq(0,2*pi,length.out = nPoints)
circleDF$lon <- unlist(lapply(centers$longitude, function(x) x + radiusLon * cos(angle)))
circleDF$lat <- unlist(lapply(centers$latitude, function(x) x + radiusLat * sin(angle)))
return(circleDF)
}
# here is the data frame for all circles
myCircles <- make_circles(data, 0.45)
##################################################################################
island = get_map(location = c(lon = -63.247593, lat = 17.631598), zoom = 13, maptype = "satellite")
islandMap = ggmap(island, extent = "panel", legend = "bottomright")
RL = geom_point(aes(x = longitude, y = latitude), data = data, color = "#ff0000")
islandMap + RL +
scale_x_continuous(limits = c(-63.280, -63.21), expand = c(0, 0)) +
scale_y_continuous(limits = c(17.605, 17.66), expand = c(0, 0)) +
########### add circles
geom_polygon(data = myCircles, aes(lon, lat, group = ID), color = "red", alpha = 0)
Ответ 2
Хорошо, поскольку ссылка на публикацию уже предлагает - переключиться на проекцию, которая основана на метрах, а затем обратно:
library(rgeos)
library(sp)
d <- SpatialPointsDataFrame(coords = data[, -1],
data = data,
proj4string = CRS("+init=epsg:4326"))
d_mrc <- spTransform(d, CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m [email protected] +no_defs"))
Теперь width
можно указать в метрах:
d_mrc_bff_mrc <- gBuffer(d_mrc, byid = TRUE, width = 450)
Измените его и добавьте в график с помощью geom_path
:
d_mrc_bff <- spTransform(d_mrc_bff_mrc, CRS("+init=epsg:4326"))
d_mrc_bff_fort <- fortify(d_mrc_bff)
islandMap +
RL +
geom_path(data=d_mrc_bff_fort, aes(long, lat, group=group), color="red") +
scale_x_continuous(limits = c(-63.280, -63.21), expand = c(0, 0)) +
scale_y_continuous(limits = c(17.605, 17.66), expand = c(0, 0))
![введите описание изображения здесь]()
Ответ 3
Расчет расстояния в км при заданной широте и долготе не является сверхпростым; Например, 1 градус лат/длиннее - это большее расстояние на экваторе, чем у полюсов. Если вы хотите легко обходное решение, которое вы можете наблюдать за точностью, вы можете попробовать:
islandMap + RL +
scale_x_continuous(limits = c(-63.280, -63.21), expand = c(0, 0)) +
scale_y_continuous(limits = c(17.605, 17.66), expand = c(0, 0)) +
geom_point(aes(x = longitude, y = latitude), data = data, size = 20, shape = 1, color = "#ff0000")
![введите описание изображения здесь]()
Вам нужно настроить параметр size
во втором geom_point
, чтобы приблизиться к тому, что вы хотите. Надеюсь, это поможет!
Ответ 4
Точное решение использует функцию geosphere:: destPoint(). Это работает без переключения проекций.
Определите функцию для определения 360 точек с определенным радиусом вокруг одной точки:
library(dplyr)
library(geosphere)
fn_circle <- function(id1, lon1, lat1, radius){
data.frame(ID = id1, degree = 1:360) %>%
rowwise() %>%
mutate(lon = destPoint(c(lon1, lat1), degree, radius)[1]) %>%
mutate(lat = destPoint(c(lon1, lat1), degree, radius)[2])
}
Применить функцию к каждой строке data
и преобразовать в data.frame:
circle <- apply(data, 1, function(x) fn_circle(x[1], x[2], x[3], 450))
circle <- do.call(rbind, circle)
Тогда отображение легко получить:
islandMap +
RL +
scale_x_continuous(limits = c(-63.280, -63.21), expand = c(0, 0)) +
scale_y_continuous(limits = c(17.605, 17.66), expand = c(0, 0)) +
geom_polygon(data = circle, aes(lon, lat, group = ID), color = "red", alpha = 0)
![введите описание изображения здесь]()