Is there a way to add a scale bar (for linear distances) to ggmap?

Question

Not that it's critical to my question, but here is my plot example, on top of which I'd like to add a scale bar.

ggmap(get_map(location = "Kinston, NC", zoom = 12, maptype = 'hybrid')) +
geom_point(x = -77.61198, y = 35.227792, colour = "red", size = 5) +
geom_point(x = -77.57306, y = 35.30288, colour = "blue", size = 3) +
geom_point(x = -77.543, y = 35.196, colour = "blue", size = 3) +
geom_text(x = -77.575, y = 35.297, label = "CRONOS Data") +
geom_text(x = -77.54, y = 35.19, label = "NOAA") +
geom_text(x = -77.61, y = 35.22, label = "PP Site")

Liverwurst

Hi, have a look at ?OSM_scale_lookup and the associated FAQ links – Felipa

The package ggsn by Oswaldo Santos adds a scale bar to maps created with ggplot or ggmap. Also, you can look to other options proposed for a similar question stackoverflow – Gstring

Answer 1

There are a few things you need to do to make this happen.

First is to put your data into a data.frame():

sites.data = data.frame(lon = c(-77.61198, -77.57306, -77.543),
                        lat = c(35.227792, 35.30288, 35.196),
                        label = c("PP Site","NOAA", "CRONOS Data"),
                        colour = c("red","blue","blue"))

Now we can get the map for this region using the gg_map package:

require(gg_map)
map.base <- get_map(location = c(lon = mean(sites.data$lon),
                                 lat = mean(sites.data$lat)),
                    zoom = 10) # could also use zoom = "auto"

We'll need the extents of that image:

bb <- attr(map.base,"bb")

Now we start figuring out the scale. First, we need a function give us the distance between two points, based on lat/long. For that, we use the Haversine formula, described by Floris at Calculate distance in (x, y) between two GPS-Points:

distHaversine <- function(long, lat){

  long <- long*pi/180
  lat <- lat*pi/180  
  dlong = (long[2] - long[1])
  dlat  = (lat[2] - lat[1])

  # Haversine formula:
  R = 6371;
  a = sin(dlat/2)*sin(dlat/2) + cos(lat[1])*cos(lat[2])*sin(dlong/2)*sin(dlong/2)
  c = 2 * atan2( sqrt(a), sqrt(1-a) )
  d = R * c
  return(d) # in km
}

The next step is to work out the points that will define our scale bar. For this example, I put something in the lower left of the plot, using the bounding box that we've already figured out:

sbar <- data.frame(lon.start = c(bb$ll.lon + 0.1*(bb$ur.lon - bb$ll.lon)),
                   lon.end = c(bb$ll.lon + 0.25*(bb$ur.lon - bb$ll.lon)),
                   lat.start = c(bb$ll.lat + 0.1*(bb$ur.lat - bb$ll.lat)),
                   lat.end = c(bb$ll.lat + 0.1*(bb$ur.lat - bb$ll.lat)))

sbar$distance = distHaversine(long = c(sbar$lon.start,sbar$lon.end),
                              lat = c(sbar$lat.start,sbar$lat.end))

Finally, we can draw the map with the scale.

ptspermm <- 2.83464567  # need this because geom_text uses mm, and themes use pts. Urgh.

map.scale <- ggmap(map.base,
                   extent = "normal", 
                   maprange = FALSE) %+% sites.data +
  geom_point(aes(x = lon,
                 y = lat,
                 colour = colour)) +
  geom_text(aes(x = lon,
                y = lat,
                label = label),
            hjust = 0,
            vjust = 0.5,
            size = 8/ptspermm) +    
  geom_segment(data = sbar,
               aes(x = lon.start,
                   xend = lon.end,
                   y = lat.start,
                   yend = lat.end)) +
  geom_text(data = sbar,
            aes(x = (lon.start + lon.end)/2,
           y = lat.start + 0.025*(bb$ur.lat - bb$ll.lat),
           label = paste(format(distance, 
                                digits = 4,
                                nsmall = 2),
                         'km')),
           hjust = 0.5,
           vjust = 0,
           size = 8/ptspermm)  +
  coord_map(projection="mercator",
            xlim=c(bb$ll.lon, bb$ur.lon),
            ylim=c(bb$ll.lat, bb$ur.lat))  

Then we save it...

# Fix presentation ----
map.out <- map.scale +  
  theme_bw(base_size = 8) +
  theme(legend.justification=c(1,1), 
        legend.position = c(1,1)) 

ggsave(filename ="map.png", 
       plot = map.out,
       dpi = 300,
       width = 4, 
       height = 3,
       units = c("in"))

Which gives you something like this:

The nice thing is that all of the plotting uses ggplot2(), so you can use the documentation at http://ggplot2.org to make this look how you need.

Enwreathe

This is an extremely helpful answer @AndyClifton! – Refrain

trying to use your code and am not sure if there is not an error somehwere? double checking your lat/long example here: [link] (movable-type.co.uk/scripts/latlong.html gives ) gives me 12.02 km instead of the 13.52 your formula gets? – Rebane

Updated November 30, 2015 to fix the problem that @maja noticed, which was caused by me not expressing the lat/Lon in radians in my distance calculation. – Enwreathe

Answer 2

I've reworked @Andy Clifton's code to add a more precise measure of the distance, and to allow for the scale bar to be of a desired length, as opposed to depending on the positioning of the bar.

Andy's code got me 99 percent of the way, but the Haversine formula used in his code is not validated with results from other sources, although I can't find the error myself.

This first part is copied from Andy Clifton's answer above just for completeness of the code:

sites.data = data.frame(lon = c(-77.61198, -77.57306, -77.543),
                        lat = c(35.227792, 35.30288, 35.196),
                        label = c("PP Site","NOAA", "CRONOS Data"),
                        colour = c("red","blue","blue"))
map.base <- get_map(location = c(lon = mean(sites.data$lon),
                                   lat = mean(sites.data$lat)),
                      zoom = 10)
bb <- attr(map.base,"bb")
sbar <- data.frame(lon.start = c(bb$ll.lon + 0.1*(bb$ur.lon - bb$ll.lon)),
                     lon.end = c(bb$ll.lon + 0.25*(bb$ur.lon - bb$ll.lon)),
                     lat.start = c(bb$ll.lat + 0.1*(bb$ur.lat - bb$ll.lat)),
                     lat.end = c(bb$ll.lat + 0.1*(bb$ur.lat - bb$ll.lat)))

The next two steps are different:

First use the distVincentyEllipsoid function from the geosphere package to calculate the distance even more preciseley than the Haversine formula:

sbar$distance <- geosphere::distVincentyEllipsoid(c(sbar$lon.start,sbar$lat.start),
 c(sbar$lon.end,sbar$lat.end))

Then correct the scale bar so that is a standard length - depending on the scale of your map. In this example 20km seems like a nice reasonable choice, i.e. 20,000 meters:

scalebar.length <- 20
sbar$lon.end <- sbar$lon.start +
((sbar$lon.end-sbar$lon.start)/sbar$distance)*scalebar.length*1000

Again using Andy's code, I've only added the arrows to the geom_segment because I think it looks nicer

ptspermm <- 2.83464567  # need this because geom_text uses mm, and themes use pts. Urgh.

map.scale <- ggmap(map.base,
                   extent = "normal", 
                   maprange = FALSE) %+% sites.data +
  geom_point(aes(x = lon,
                 y = lat,
                 colour = colour)) +
  geom_text(aes(x = lon,
                y = lat,
                label = label),
            hjust = 0,
            vjust = 0.5,
            size = 8/ptspermm) +    
  geom_segment(data = sbar,
               aes(x = lon.start,
                   xend = lon.end,
                   y = lat.start,
                   yend = lat.end),
               arrow=arrow(angle = 90, length = unit(0.1, "cm"),
                           ends = "both", type = "open")) +
  geom_text(data = sbar,
            aes(x = (lon.start + lon.end)/2,
                y = lat.start + 0.025*(bb$ur.lat - bb$ll.lat),
                label = paste(format(scalebar.length),
                              'km')),
            hjust = 0.5,
            vjust = 0,
            size = 8/ptspermm)  +
  coord_map(projection = "mercator",
            xlim=c(bb$ll.lon, bb$ur.lon),
            ylim=c(bb$ll.lat, bb$ur.lat))  

# Fix presentation ----
map.out <- map.scale +  
  theme_bw(base_size = 8) +
  theme(legend.justification = c(1,1), 
        legend.position = c(1,1)) 

ggsave(filename ="map.png", 
       plot = map.out,
       dpi = 300,
       width = 4, 
       height = 3,
       units = c("in"))

Rebane