Нарисуйте круг с определенным радиусом на карте - убедитесь, что расстояние правильное

В R я пытаюсь нарисовать круг на карте с определенным расстоянием, учитывая его центроид как широту, длинные точки. Я нашел эту тему и использовал код @lukeA для достижения того, чего я хочу. Однако, кажется, что расстояние не правильное. Расстояние, которое я получаю между двумя долготами на определенной широте, не соответствует тому, что изображено на графике. Веб-сайт, который можно использовать для измерения расстояния: http://www.stevemorse.org/nearest/distance.php

Ниже приведен код. Единицы в spTransform указаны как us-mi, поэтому мы должны указать диаметр окружности в милях.

Мне нужен круг с центром в широте = 45, долготе = -90. Расстояние до широты = 45, долготы = -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") 

Карта с кругом


person ilyas    schedule 19.04.2017    source источник
comment
Я не думаю, что вам нужно умножать радиус на 2. Поскольку в документации gBuffer указано, что ширина - это расстояние от исходной геометрии для включения в новую геометрию.   -  person ahly    schedule 19.04.2017
comment
Тогда это слишком мало.   -  person ilyas    schedule 19.04.2017
comment
Попробуйте изменить CRS следующим образом: d_mrc <- spTransform(d, CRS("+proj=utm +zone=16 +datum=WGS84 +units=us-mi")).   -  person ahly    schedule 19.04.2017
comment
@ahly Кажется, это работает. Как вы нашли эту информацию? Где я могу проверить правильную зону, проекцию и т.д.?   -  person ilyas    schedule 19.04.2017
comment
Вы можете рассчитать зону, используя эту функцию: function(long) { (floor((long + 180)/6) %% 60) + 1 }. Для получения дополнительной информации о проекционных системах вы можете обратиться по адресу nceas.ucsb.edu/~ Фрейзер/RSpatialGuides/   -  person ahly    schedule 19.04.2017
comment
Возможно, взгляните на geosphere::destPoint(), как описано (например) в vignette("geosphere", package="geosphere")   -  person Josh O'Brien    schedule 19.04.2017
comment
@ahly, можете ли вы представить свое решение в качестве ответа, чтобы я мог его принять. Спасибо.   -  person ilyas    schedule 21.04.2017
comment
@ilyas добавил ответ.   -  person ahly    schedule 22.04.2017


Ответы (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") 

введите здесь описание изображения

person ahly    schedule 21.04.2017