在地图上绘制具有一定半径的圆 - 确保距离正确
Draw circle with certain radius on a map - Making sure distance is right
在 R 中,我试图在地图上绘制一个具有一定距离的圆,给定其质心作为纬度、经度点。我发现 and used @lukeA
's code to achieve what I want. However, it seems like the distance is not right. The distance I get between two longitudes at a certain latitude does not correspond to what is plotted. The website that one can use to measure distance is: http://www.stevemorse.org/nearest/distance.php
下面是代码。 spTransform
中的单位指定为 us-mi
,因此我们应该能够给出以英里为单位的圆直径。
我想要一个以纬度 = 45,经度 = -90 为中心的圆。到 lat = 45, long = -86 的距离约为 195 英里。因此,直径等于 2*195=390 英里。但是画出的圆圈很远。问题可能与投影有关。
有人可以帮助我解决我所缺少的问题吗?
library(tidyverse)
library(data.table)
library(rgeos)
library(sp)
states = c('illinois', 'wisconsin', 'iowa', 'minnesota')
state_4s = map_data('state') %>% data.table() %>% subset(region %in% states)
counties_4s = counties[region %in% states ]
map_4s <- ggplot(data = state_4s,
mapping = aes(x = long, y = lat, group = group)) +
geom_polygon(color = "black", fill = "gray", size = 0.1) +
theme_bw()
map_4s = map_4s +
geom_polygon(data = counties_4s, fill = NA, color = 'white', size = 0.1) +
geom_polygon(color = "black", fill = NA, size = 0.1) +
coord_map("mercator")
long = c(-90)
lat = c(45)
center = data.frame(long, lat)
d <- SpatialPointsDataFrame(coords = center,
data = center,
proj4string = CRS("+init=epsg:4326"))
d_mrc <- spTransform(d, CRS("+proj=merc +init=epsg:4326 +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k_0=1.0 +units=us-mi +nadgrids=@null +no_defs"))
# Now, the width can be specified in miles:
d_mrc_bff_mrc <- gBuffer(d_mrc, byid = T, width = 195*2, capStyle = 'round')
d_mrc_bff <- spTransform(d_mrc_bff_mrc, CRS("+init=epsg:4326"))
d_mrc_bff_fort <- fortify(d_mrc_bff)
map_4s +
geom_point(data = center, aes(x = long, y = lat, group = 1)) +
geom_path(data=d_mrc_bff_fort, aes(long, lat, group=group), color="red")
library(tidyverse)
library(data.table)
library(rgeos)
library(sp)
states = c('illinois', 'wisconsin', 'iowa', 'minnesota')
state_4s = map_data('state') %>% data.table() %>% subset(region %in% states)
map_4s <- ggplot(data = state_4s,
mapping = aes(x = long, y = lat, group = group)) +
geom_polygon(color = "black", fill = "gray", size = 0.1) +
theme_bw()
long = c(-90)
lat = c(45)
center = data.frame(long, lat)
d <- SpatialPointsDataFrame(coords = center,
data = center,
proj4string = CRS("+init=epsg:4326"))
long2UTMZone <- function(long) { (floor((long + 180)/6) %% 60) + 1 }
d_mrc <- spTransform(d, CRS(paste0("+proj=utm +zone=", long2UTMZone(long)," +datum=WGS84 +units=us-mi")))
# Now, the width can be specified in miles:
d_mrc_bff_mrc <- gBuffer(d_mrc, byid = T, width = 195, capStyle = 'round')
d_mrc_bff <- spTransform(d_mrc_bff_mrc, CRS("+init=epsg:4326"))
d_mrc_bff_fort <- fortify(d_mrc_bff)
map_4s +
geom_point(data = center, aes(x = long, y = lat, group = 1)) +
geom_path(data=d_mrc_bff_fort, aes(long, lat, group=group), color="red")
在 R 中,我试图在地图上绘制一个具有一定距离的圆,给定其质心作为纬度、经度点。我发现 @lukeA
's code to achieve what I want. However, it seems like the distance is not right. The distance I get between two longitudes at a certain latitude does not correspond to what is plotted. The website that one can use to measure distance is: http://www.stevemorse.org/nearest/distance.php
下面是代码。 spTransform
中的单位指定为 us-mi
,因此我们应该能够给出以英里为单位的圆直径。
我想要一个以纬度 = 45,经度 = -90 为中心的圆。到 lat = 45, long = -86 的距离约为 195 英里。因此,直径等于 2*195=390 英里。但是画出的圆圈很远。问题可能与投影有关。
有人可以帮助我解决我所缺少的问题吗?
library(tidyverse)
library(data.table)
library(rgeos)
library(sp)
states = c('illinois', 'wisconsin', 'iowa', 'minnesota')
state_4s = map_data('state') %>% data.table() %>% subset(region %in% states)
counties_4s = counties[region %in% states ]
map_4s <- ggplot(data = state_4s,
mapping = aes(x = long, y = lat, group = group)) +
geom_polygon(color = "black", fill = "gray", size = 0.1) +
theme_bw()
map_4s = map_4s +
geom_polygon(data = counties_4s, fill = NA, color = 'white', size = 0.1) +
geom_polygon(color = "black", fill = NA, size = 0.1) +
coord_map("mercator")
long = c(-90)
lat = c(45)
center = data.frame(long, lat)
d <- SpatialPointsDataFrame(coords = center,
data = center,
proj4string = CRS("+init=epsg:4326"))
d_mrc <- spTransform(d, CRS("+proj=merc +init=epsg:4326 +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k_0=1.0 +units=us-mi +nadgrids=@null +no_defs"))
# Now, the width can be specified in miles:
d_mrc_bff_mrc <- gBuffer(d_mrc, byid = T, width = 195*2, capStyle = 'round')
d_mrc_bff <- spTransform(d_mrc_bff_mrc, CRS("+init=epsg:4326"))
d_mrc_bff_fort <- fortify(d_mrc_bff)
map_4s +
geom_point(data = center, aes(x = long, y = lat, group = 1)) +
geom_path(data=d_mrc_bff_fort, aes(long, lat, group=group), color="red")
library(tidyverse)
library(data.table)
library(rgeos)
library(sp)
states = c('illinois', 'wisconsin', 'iowa', 'minnesota')
state_4s = map_data('state') %>% data.table() %>% subset(region %in% states)
map_4s <- ggplot(data = state_4s,
mapping = aes(x = long, y = lat, group = group)) +
geom_polygon(color = "black", fill = "gray", size = 0.1) +
theme_bw()
long = c(-90)
lat = c(45)
center = data.frame(long, lat)
d <- SpatialPointsDataFrame(coords = center,
data = center,
proj4string = CRS("+init=epsg:4326"))
long2UTMZone <- function(long) { (floor((long + 180)/6) %% 60) + 1 }
d_mrc <- spTransform(d, CRS(paste0("+proj=utm +zone=", long2UTMZone(long)," +datum=WGS84 +units=us-mi")))
# Now, the width can be specified in miles:
d_mrc_bff_mrc <- gBuffer(d_mrc, byid = T, width = 195, capStyle = 'round')
d_mrc_bff <- spTransform(d_mrc_bff_mrc, CRS("+init=epsg:4326"))
d_mrc_bff_fort <- fortify(d_mrc_bff)
map_4s +
geom_point(data = center, aes(x = long, y = lat, group = 1)) +
geom_path(data=d_mrc_bff_fort, aes(long, lat, group=group), color="red")