Czy istnieje sposób dodania paska skali (dla odległości liniowych) do ggmap?

Question

Nie ma to krytycznego znaczenia dla mojego pytania, ale tutaj jest mój przykład fabuły, nad którym chciałbym dodać pasek skali.

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

Istnieje kilka rzeczy, które trzeba zrobić, aby tak się stało.

Pierwszy jest umieścić swoje dane w 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")) 

Teraz możemy dostać mapę tego regionu, korzystając z pakietu 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" 

Będziemy potrzebować stopnie tego image:

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

Teraz zaczynamy obliczać skalę. Po pierwsze, potrzebujemy funkcji, która daje nam odległość między dwoma punktami, na podstawie lat/długo. W tym celu używamy formuły Haversine, opisanej przez Floris w 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 
} 

Następnym krokiem jest wypracowanie punktów, które definiują nasz pasek skali. W tym przykładzie, umieścić coś w lewym dolnym rogu wykresu, za pomocą obwiedni, że mamy już zorientowali się:

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

Wreszcie możemy narysować mapę z skali.

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

Potem go uratować ...

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

co daje mniej więcej tak:

Dobrą rzeczą jest to, że cały spisek używa ggplot2(), więc możesz skorzystać z dokumentacji pod numerem http://ggplot2.org, aby wyglądała tak, jak potrzebujesz.

Answer 2

Przerobiłem kod @Andy Clifton, aby dodać dokładniejszy pomiar odległości i aby pasek skali miał pożądaną długość, w przeciwieństwie do zależności od położenia paska.

Kod Andy'ego dostarczył mi 99 procent drogi, ale formuła Haversine użyta w jego kodzie nie jest sprawdzana z wynikami z innych źródeł, chociaż sam nie mogę znaleźć błędu.

Ta pierwsza część jest kopiowana z odpowiedzią Andy Cliftona powyżej tylko dla kompletności kodu:

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

Kolejne dwa kroki są różne:

Najpierw użyj funkcji distVincentyEllipsoid z pakietu geosphere obliczyć odległość jeszcze bardziej precyzyjnie niż formuła Haversine:

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

Następnie popraw podziałkę skali więc jest to standardowa długość - w zależności od skali mapy. W tym przykładzie 20 km wydaje się miłym rozsądnego wyboru, czyli 20.000 metrów:

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

ponownie stosując kod Andy'ego, Dodałem tylko strzałek do geom_segment bo myślę, że wygląda ładniej

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