mexico - maps in r ggplot2

¿Hay alguna manera de agregar una barra de escala(para distancias lineales) a ggmap?

Question

No es que sea crítico para mi pregunta, pero aquí está mi ejemplo de trama, sobre la cual me gustaría agregar una barra de escala.

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

Answer 1

Hay algunas cosas que debe hacer para que esto suceda.

Lo primero es poner tus datos en un 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"))

Ahora podemos obtener el mapa para esta región usando el paquete gg_map :

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"

Necesitaremos las extensiones de esa imagen:

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

Ahora empezamos a averiguar la escala. Primero, necesitamos una función que nos dé la distancia entre dos puntos, en función de lat / long. Para eso, utilizamos la fórmula de Haversine, descrita por Floris en la distancia calculada en (x, y) entre dos puntos GPS :

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
}

El siguiente paso es resolver los puntos que definirán nuestra barra de escala. Para este ejemplo, coloco algo en la parte inferior izquierda de la gráfica, utilizando el cuadro delimitador que ya hemos descubierto:

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

Finalmente, podemos dibujar el mapa con la escala.

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

Luego lo guardamos ...

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

Lo que te da algo como esto:

Lo bueno es que todo el trazado utiliza ggplot2() , por lo que puede usar la documentación en http://ggplot2.org para hacer que se vea como lo necesita.

Answer 2

He modificado el código de @Andy Clifton para agregar una medida más precisa de la distancia y permitir que la barra de escala tenga la longitud deseada, en lugar de depender de la posición de la barra.

El código de Andy me dio el 99% del camino, pero la fórmula de Haversine utilizada en su código no está validada con los resultados de otras fuentes, aunque yo mismo no puedo encontrar el error.

Esta primera parte se copia de la respuesta de Andy Clifton que aparece arriba para completar el código:

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

Los siguientes dos pasos son diferentes:

Primero use la función distVincentyEllipsoid del paquete de geosphere para calcular la distancia aún más precisa que la fórmula de Haversine:

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

Luego corrija la barra de escala para que sea una longitud estándar, dependiendo de la escala de su mapa. En este ejemplo, 20 km parece una buena opción razonable, es decir, 20,000 metros:

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

De nuevo, usando el código de Andy, solo agregué las flechas a geom_segment porque creo que se ve mejor

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