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")
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.
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"))
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With