Plot circle with a certain radius around point on a map in ggplot2

Question

I have a map with the 8 points plotted on it:

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

Now I want to plot a circle around each of the 8 plotted locations. The circle has to have a radius of 450 meters.

This is what I mean, but then using ggplot: https://gis.stackexchange.com/questions/119736/ggmap-create-circle-symbol-where-radius-represents-distance-miles-or-km

How can I achieve this?

Hormone

Answer 1

If you only work on a small area of the earth, here is a approximation. Each degree of the latitude represents 40075 / 360 kilometers. Each degrees of longitude represents (40075 / 360) * cos(latitude) kilomemters. With this, we can calculate approximately a data frame including all points on circles, knowing the circle centers and radius.

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)

Unvoice

Really nicely done. Follow the 111, but where does the 57.3 come in? – Tempered

1 radian = 57.3 degree (180 / pi) – Unvoice

Answer 2

Well, as the referred posting already suggests - switch to a projection that is based in meters, and then back:

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 +nadgrids=@null +no_defs"))

Now, the width can be specified in meters:

d_mrc_bff_mrc <- gBuffer(d_mrc, byid = TRUE, width = 450)

Transform it back and add it to the plot using 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)) 

Preordain

Answer 3

Calculating distance in km given latitude and longitude isn't super straightforward; 1 degree lat/long is a greater distance at the equator than at the poles, for example. If you want an easy workaround that you can eyeball for accuracy, you might try:

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

You'll need to adjust the size paramter in the 2nd geom_point to get closer to what you want. I hope that helps!

Wulfila

Thank you very much! This could be a solution, however it is not as accurate as I hoped for.. Any other suggestions are more than welcome. – Hormone

Definitely; it isn't visually reliable, and shouldn't be treated as such if you publish for example. I'm not sure what the relationship is between the size of a point and the axis marks. I've tried to do something similar though, and because the Earth isn't flat (apparently) you have to incorporate a z-axis in order to accurately show a circle on a flat map. – Wulfila

Also note that if you make large enough "circles" they'll get stretched into ovals, once again because of sphere stuff. It's a bummer, really haha. – Wulfila

Thanks Nancy! Would have been so much easier if we stayed with 'the earth is flat'-thing.. – Hormone

Answer 4

An accurate solution is using the geosphere::destPoint() function. This works without switching projections.

Define function to determine 360 points with a certain radius around one point:

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]) 
}

Apply function to each row of data and convert to data.frame:

circle <- apply(data, 1, function(x) fn_circle(x[1], x[2], x[3], 450))
circle <- do.call(rbind, circle)

Then the map can be easily obtained by:

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)

Fastigiate

Answer 5

A solution using st_buffer() from the sf package.

library(ggmap)
library(ggplot2)
library(sf)

data <- data.frame(
  ID = 1:8,
  longitude = c(-63.27462, -63.26499, -63.25658, -63.2519, 
                -63.2311, -63.2175, -63.23623, -63.25958),
  latitude = c(17.6328, 17.64614, 17.64755, 17.64632, 
               17.64888, 17.63113, 17.61252, 17.62463)
)

Convert data.frame to sf object:

points_sf <- sf::st_as_sf(data, coords = c("longitude", "latitude"), crs = 4326)

For this example we use UTM zone 20, which contains the coordinates of the island:

data_sf_utm <- sf::st_transform(points_sf, "+proj=utm +zone=20")

Now we can buffer the point by 450 meters:

circle <- sf::st_buffer(data_sf_utm, dist = 450)

ggmap seems to have some issues with geom_sf. Setting inherit.aes to FALSE returns the desired map.

island <- ggmap::get_map(location = c(lon = -63.247593, lat = 17.631598), zoom = 14, maptype = "satellite")

ggmap(island, extent = "panel", legend = "bottomright") + 
  geom_sf(data = points_sf, color = "red", inherit.aes = FALSE) +
  geom_sf(data = circle, color = "red", alpha = 0, inherit.aes = FALSE)

^{Created on 2020-10-11 by the reprex package (v0.3.0)}

Fastigiate