Draw circle with certain radius on a map - Making sure distance is right

Question

In R, I am trying to plot a circle on a map with certain distance given its centroid as lat, long points. I found this thread 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

Below is the code. The units in spTransform is specified as us-mi, so we should be able to give circle diameter in miles.

I want a circle centered at lat = 45, long = -90. The distance to lat = 45, long = -86 is ~195 mi. The diameter thus equals to 2*195=390 mi. But the circle drawn goes far out. The problem might be related to projections.

Can someone help me what I'm missing?

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

I don't think you have to multiply the radius by 2. As the documentation of `gBuffer` states width is "Distance from original geometry to include in the new geometry". — ahly, Apr 19 '17 at 19:43
Try changing your CRS like this: `d_mrc <- spTransform(d, CRS("+proj=utm +zone=16 +datum=WGS84 +units=us-mi"))`. — ahly, Apr 19 '17 at 19:56
@ahly This seems to work. How did you find this info? Where can I check out proper zone, projection etc.? — ilyas, Apr 19 '17 at 20:04
You can calculate the zone using this function: `function(long) { (floor((long + 180)/6) %% 60) + 1 }`. For more information about projection systems, you can refer https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf — ahly, Apr 19 '17 at 20:22
Perhaps have a look at `geosphere::destPoint()`, as documented (for example) in `vignette("geosphere", package="geosphere")` — Josh O'Brien, Apr 19 '17 at 20:33
@ahly, can you submit your solution as an answer, so I can accept it. Thanks. — ilyas, Apr 21 '17 at 16:42

score 0 · Accepted Answer · answered Apr 21 '17 at 19:07

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

Draw circle with certain radius on a map - Making sure distance is right

1 Answers1