在地图上绘制具有一定半径的圆 - 确保距离是正确的答案

【问题标题】：Draw circle with certain radius on a map - Making sure distance is right在地图上绘制具有一定半径的圆 - 确保距离是正确的
【发布时间】：2017-04-19 19:16:42
【问题描述】：

在 R 中，我试图在地图上绘制一个具有一定距离的圆，因为它的质心为 lat，long 点。我找到了this thread 并使用@lukeA 的代码来实现我想要的。但是，距离似乎不对。我在某个纬度的两个经度之间得到的距离与绘制的不对应。可以用来测量距离的网站是：http://www.stevemorse.org/nearest/distance.php

下面是代码。 spTransform 中的单位指定为us-mi，因此我们应该能够以英里为单位给出圆直径。

我想要一个以 lat = 45、long = -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")

【问题讨论】：

我认为您不必将半径乘以 2。正如 gBuffer 的文档所述，宽度是“与原始几何图形的距离以包含在新几何图形中”。
那就太小了。
尝试像这样更改您的 CRS：d_mrc <- spTransform(d, CRS("+proj=utm +zone=16 +datum=WGS84 +units=us-mi")).
@ahly 这似乎有效。你是怎么找到这个信息的？我在哪里可以查看正确的区域、投影等？
您可以使用此函数计算区域：function(long) { (floor((long + 180)/6) %% 60) + 1 }。有关投影系统的更多信息，您可以参考nceas.ucsb.edu/~frazier/RSpatialGuides/…

标签： r ggplot2 map-projections

【解决方案1】：

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

【讨论】：