Hot questions for Using Ggplot2 in ggmap

Question:

I want to get a map with RgoogleMaps from R, with a specific coordinates boundary.

What I can call is GetMap, and specify a center, I must add a zoom level. Everything works fine, except that I am not getting an image map bounded with the coordinates I choose.

Here's an example:

lat <- c(44.49,44.5)                
lon <- c(11.33,11.36)               
center = c(mean(lat), mean(lon))    
zoom <- 14                          
mmap <- GetMap(center = center, zoom=zoom, maptype= "satellite", destfile = "m.png") 

The problem is that only the center is passed as a parameter, and thus the whole image I see is dependant on the zoom level. So, I cannot really understand what are the boundaries of the image I get. What I want to do is to get an image bounded exactly with the coordinates I am defining. Is this possible (also with other map packages)?


Answer:

Here is one way. First, you get a map with a certain zoom. Then, you add the lon and lat limit when you draw a figure, which you can do with scale_x_continuous and scale_y_continuous.

library(ggmap)
library(ggplot2)

### Set a range
lat <- c(44.49, 44.5)                
lon <- c(11.33, 11.36)   

### Get a map
map <- get_map(location = c(lon = mean(lon), lat = mean(lat)), zoom = 14,
               maptype = "satellite", source = "google")

### When you draw a figure, you limit lon and lat.      
foo <- ggmap(map)+
       scale_x_continuous(limits = c(11.33, 11.36), expand = c(0, 0)) +
       scale_y_continuous(limits = c(44.49, 44.5), expand = c(0, 0))

foo

Question:

I am generating maps with world-scale data, and then zooming in to certain regions. On the zoomed-in view, I would like to show that there are other data points outside the bounding box, by putting arrowheads that point from the center of the box to where the data point is in the outside world.

Note: I do not need it to be a "great circle" path, just XY vectors in Mercator projection, because I imagine this will be useful for "normal" plots as well.

As an example, here is the world map showing the extent of the data:

And here is the zoomed in view, with magenta arrows manually added to show what I would like to generate.

Below is the code and data I am using to generate these two basic plots. What I need is a way to generate the arrowheads.

require(ggplot2)

te = structure(list(lat = c(33.7399, 32.8571, 50.2214, 36.96263, 33.5835, 
33.54557, 47.76147, 48, 59.40289, 35.93411, 32.87962, 38.3241, 
50.03844, 37.44, 50.07774, 50.26668, 36.5944), lng = c(-118.37608, 
-117.25746, -5.3865, -122.00809, -117.86159, -117.79805, -124.45055, 
-126, -146.35157, -122.931472, -117.25285, -123.07331, -5.26339, 
25.4, -5.709894, -3.86828, -121.96201)), .Names = c("lat", "lng"
), class = "data.frame", row.names = c(NA, -17L))

all_states = map_data("world")

# world version:
wp = ggplot() + 
      geom_polygon(data = all_states, aes(x = long, y = lat, group = group), colour = "gray",
                   fill = "gray") +
      coord_cartesian(ylim = c(0, 80), xlim = c(-155, 45)) + 
      geom_point(data = te, aes(x = lng, y = lat), color = "blue", size = 5,alpha = 0.6)

print(wp)

#states plot
sp = ggplot() +
      geom_polygon(data = all_states, aes(x = long, y = lat, group = group), colour = "gray", fill = "gray") +
      coord_cartesian(ylim = c(30, 52), xlim = c(-128, -114)) + 
      geom_point(data = te, aes(x = lng, y = lat), color = "blue", size = 5, alpha = 0.6) 

print(sp)

Answer:

This solution uses sp and rgeos packages to manipulate spatial data, the main crux being intersecting lines and a box polygon to get the edge points for arrows. Then if you draw arrows with geom_segment and zero width, the line is invisible and only the arrow head remains.

This function computes the line-box intersections:

boxint <- function(xlim, ylim, xp, yp){
    ## build box as SpatialPolygons
    box = cbind(xlim[c(1,2,2,1,1)],
        ylim[c(1,1,2,2,1)])
    box <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(box)),ID=1)))

    ## get centre of box
    x0=mean(xlim)
    y0=mean(ylim)

    ## construct line segments to points
    sl = sp::SpatialLines(
        lapply(1:length(xp),
               function(i){
                   sp::Lines(list(sp::Line(cbind(c(x0,xp[i]),c(y0,yp[i])))),ID=i)
               }
               )
        )
    ## intersect lines segments with boxes to make points
    pts = rgeos::gIntersection(sl, as(box, "SpatialLines"))
    as.data.frame(sp::coordinates(pts), row.names=1:length(xp))
}

And this returns the geom with arrows:

wherelse <- function(xlim, ylim, points){
    ## get points outside bounding box
    outsides = points[!(
        points$lng>=xlim[1] &
            points$lng <= xlim[2] &
                points$lat >= ylim[1] &
                    points$lat <= ylim[2]),]
    npts = nrow(outsides)
    ## get centre point of box
    x = rep(mean(xlim),npts)
    y = rep(mean(ylim),npts)

    ## compute box-point intersections
    pts = boxint(xlim, ylim, outsides$lng, outsides$lat)
    pts$x0=x
    pts$y0=y
    ## create arrow segments as invisible lines with visible arrowheads
    ggplot2::geom_segment(data=pts, aes(x=x0,y=y0,xend=x,yend=y),
       lwd=0, arrow=grid::arrow(length=unit(0.5,"cm"),
       type="closed"),col="magenta")
}

So your example, the basic plot is:

sp = ggplot() + 
  geom_polygon(
   data=all_states, 
    aes(x=long, y=lat, group = group),colour="gray",fill="gray" ) + 
    coord_cartesian(ylim=c(30, 52), xlim=c(-128,-114)) + 
    geom_point(data=te,aes(x=lng,y=lat),color="blue",size=5,alpha=0.6)

and then add the arrows with:

sp + wherelse(c(-128,-114), c(30,52), te)

Not sure if there's an option to draw arrows exactly like you want them though!

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?


Answer:

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)

Question:

I am attempting to use ggmap to look at education scores by school. I created a coordinate list of all the schools and the individual student scores like so:

     score      lat       lon
3205    45 28.04096 -82.54980
8275    60 27.32163 -80.37673
4645    38 27.45734 -82.52599
8962    98 26.54113 -81.84399
9199    98 27.88948 -82.31770
340     53 26.36528 -81.79639

I first used the pattern from most of the tutorials that I worked through: http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf http://www.geo.ut.ee/aasa/LOOM02331/heatmap_in_R.html

library(ggmap)
library(RColorBrewer)

MyMap <- get_map(location = "Orlando, FL", 
                 source = "google", maptype = "roadmap", crop = FALSE, zoom = 7)

YlOrBr <- c("#FFFFD4", "#FED98E", "#FE9929", "#D95F0E", "#993404")

ggmap(MyMap) +
stat_density2d(data = s_rit, aes(x = lon, y = lat, fill = ..level.., alpha = ..level..),
               geom = "polygon", size = 0.01, bins = 16) +
scale_fill_gradient(low = "red", high = "green") +
scale_alpha(range = c(0, 0.3), guide = FALSE)

In the first plot, the graphics look great but it doesn't take the score into account.

In order to incorporate the score variable, I used this example Density2d Plot using another variable for the fill (similar to geom_tile)?:

ggmap(MyMap) %+% s_rit +
  aes(x = lon, y = lat, z = score) +
  stat_summary2d(fun = median, binwidth = c(.5, .5), alpha = 0.5) +
  scale_fill_gradientn(name = "Median", colours = YlOrBr, space = "Lab") +
  labs(x = "Longitude", y = "Latitude") +
  coord_map()

It colours by value, but it doesn't have the look of the first. The square boxes are clunky and arbitrary. Adjusting the size of the box does not help. The dispersion of the first heatmap is preferred. Is there a way to blend the look of the first graph with the value-based plot of the second?

Data

s_rit <- structure(list(score = c(45, 60, 38, 98, 98, 53, 90, 42, 96, 
45, 89, 18, 66, 2, 45, 98, 6, 83, 63, 86, 63, 81, 70, 8, 78, 
15, 7, 86, 15, 63, 55, 13, 83, 76, 78, 70, 64, 88, 61, 78, 4, 
7, 1, 70, 88, 58, 70, 58, 11, 45, 28, 42, 45, 73, 85, 86, 25, 
17, 53, 95, 49, 80, 70, 35, 94, 61, 39, 76, 28, 1, 18, 93, 73, 
67, 56, 38, 45, 66, 18, 76, 91, 76, 52, 60, 2, 38, 73, 95, 1, 
76, 6, 25, 76, 81, 35, 49, 85, 55, 66, 90), lat = c(28.040961, 
27.321633, 27.457342, 26.541129, 27.889476, 26.365284, 28.555024, 
26.541129, 26.272728, 28.279994, 27.889476, 28.279994, 26.6674, 
26.272728, 25.776045, 26.541129, 30.247658, 26.365284, 25.450123, 
27.889476, 26.541129, 27.264513, 26.718652, 28.044369, 28.251435, 
27.264513, 26.272728, 26.272728, 28.040961, 30.312239, 27.889476, 
26.541129, 26.6674, 27.321633, 26.365284, 28.279994, 26.718652, 
30.23286, 28.040961, 30.193704, 30.312239, 28.044369, 27.457342, 
25.450123, 30.23286, 30.312239, 30.193704, 28.279994, 30.247658, 
26.541129, 26.365284, 28.279994, 27.321633, 25.776045, 26.272728, 
30.23286, 30.312239, 26.718652, 26.541129, 25.450123, 28.251435, 
28.185751, 25.450123, 28.040961, 27.321633, 28.279994, 27.321633, 
27.321633, 27.321633, 28.279994, 26.718652, 28.362308, 27.264513, 
26.365284, 28.279994, 30.23286, 25.450123, 28.362308, 25.450123, 
25.776045, 30.193704, 28.251435, 27.457342, 27.321633, 28.185751, 
27.457342, 27.889476, 26.541129, 26.541129, 30.23286, 30.312239, 
26.718652, 25.450123, 26.139258, 28.040961, 30.23286, 26.718652, 
28.185751, 28.044369, 28.555024), lon = c(-82.5498, -80.376729, 
-82.525985, -81.843986, -82.317701, -81.796389, -81.276464, -81.843986, 
-80.207508, -81.331178, -82.317701, -81.331178, -80.072089, -80.207508, 
-80.199437, -81.843986, -81.808664, -81.796389, -80.433557, -82.317701, 
-81.843986, -80.432125, -80.091078, -82.394639, -81.490407, -80.432125, 
-80.207508, -80.207508, -82.5498, -81.575916, -82.317701, -81.843986, 
-80.072089, -80.376729, -81.796389, -81.331178, -80.091078, -81.585975, 
-82.5498, -81.579846, -81.575916, -82.394639, -82.525985, -80.433557, 
-81.585975, -81.575916, -81.579846, -81.331178, -81.808664, -81.843986, 
-81.796389, -81.331178, -80.376729, -80.199437, -80.207508, -81.585975, 
-81.575916, -80.091078, -81.843986, -80.433557, -81.490407, -81.289394, 
-80.433557, -82.5498, -80.376729, -81.331178, -80.376729, -80.376729, 
-80.376729, -81.331178, -80.091078, -81.428494, -80.432125, -81.796389, 
-81.331178, -81.585975, -80.433557, -81.428494, -80.433557, -80.199437, 
-81.579846, -81.490407, -82.525985, -80.376729, -81.289394, -82.525985, 
-82.317701, -81.843986, -81.843986, -81.585975, -81.575916, -80.091078, 
-80.433557, -80.238901, -82.5498, -81.585975, -80.091078, -81.289394, 
-82.394639, -81.276464)), .Names = c("score", "lat", "lon"), row.names = c(3205L, 
8275L, 4645L, 8962L, 9199L, 340L, 5381L, 8998L, 5476L, 4956L, 
9256L, 4940L, 6681L, 5586L, 1046L, 9017L, 1878L, 323L, 4175L, 
9236L, 8968L, 6885L, 5874L, 9412L, 6434L, 7168L, 5420L, 5680L, 
3202L, 1486L, 9255L, 9009L, 6833L, 8271L, 261L, 5024L, 8028L, 
1774L, 3329L, 1824L, 1464L, 9468L, 4643L, 4389L, 1506L, 1441L, 
1826L, 4968L, 1952L, 8803L, 339L, 4868L, 8266L, 1334L, 5483L, 
1727L, 1389L, 7944L, 8943L, 4416L, 6440L, 526L, 4478L, 3117L, 
8308L, 4891L, 8290L, 8299L, 8233L, 4848L, 7922L, 5795L, 6971L, 
179L, 4990L, 1776L, 4431L, 5718L, 4268L, 1157L, 1854L, 6433L, 
4637L, 8194L, 560L, 4694L, 9274L, 8903L, 8877L, 1586L, 1398L, 
5865L, 4209L, 6075L, 3307L, 1634L, 8108L, 514L, 9453L, 5210L), class = "data.frame")

Answer:

I'd like to suggest an alternate way of visualizing the distribution of scores (in general) and the median outcomes of each school. It might be better (I don't really know your data or overall problem statement) to show the distribution of scores themselves by various levels (0-10, 10-20, etc) separately then show a view of the actual median rankings per school. Something like this:

library(ggplot2)
library(ggthemes)
library(viridis) # devtools::install_github("sjmgarnier/viridis)
library(ggmap)
library(scales)
library(grid)
library(dplyr)
library(gridExtra)

dat$cut <- cut(dat$score, breaks=seq(0,100,11), labels=sprintf("Score %d-%d",seq(0, 80, 10), seq(10,90,10)))

orlando <- get_map(location="orlando, fl", source="osm", color="bw", crop=FALSE, zoom=7)

gg <- ggmap(orlando)
gg <- gg + stat_density2d(data=dat, aes(x=lon, y=lat, fill=..level.., alpha=..level..),
                          geom="polygon", size=0.01, bins=5)
gg <- gg + scale_fill_viridis()
gg <- gg + scale_alpha(range=c(0.2, 0.4), guide=FALSE)
gg <- gg + coord_map()
gg <- gg + facet_wrap(~cut, ncol=3)
gg <- gg + labs(x=NULL, y=NULL, title="Score Distribution Across All Schools\n")
gg <- gg + theme_map(base_family="Helvetica")
gg <- gg + theme(plot.title=element_text(face="bold", hjust=1))
gg <- gg + theme(panel.margin.x=unit(1, "cm"))
gg <- gg + theme(panel.margin.y=unit(1, "cm"))
gg <- gg + theme(legend.position="right")
gg <- gg + theme(strip.background=element_rect(fill="white", color="white"))
gg <- gg + theme(strip.text=element_text(face="bold", hjust=0))
gg

median_scores <- summarise(group_by(dat, lon, lat), median=median(score))
median_scores$school <- sprintf("School #%d", 1:nrow(median_scores))

gg <- ggplot(median_scores)
gg <- gg + geom_segment(aes(x=reorder(school, median), 
                            xend=reorder(school, median), 
                            y=0, yend=median), size=0.5)
gg <- gg + geom_point(aes(x=reorder(school, median), y=median))
gg <- gg + geom_text(aes(x=reorder(school, median), y=median, label=median), size=3, hjust=-0.75)
gg <- gg + scale_y_continuous(expand=c(0, 0), limits=c(0, 100))
gg <- gg + labs(x=NULL, y=NULL, title="Median Score Per School")
gg <- gg + coord_flip()
gg <- gg + theme_tufte(base_family="Helvetica")
gg <- gg + theme(axis.ticks.x=element_blank())
gg <- gg + theme(axis.text.x=element_blank())
gg <- gg + theme(plot.title=element_text(face="bold", hjust=1))
gg_med <- gg

# tweak hjust and potentially y as needed
median_scores$hjust <- 0
median_scores[median_scores$school=="School #10",]$hjust <- 1.5
median_scores[median_scores$school=="School #8",]$hjust <- 0
median_scores[median_scores$school=="School #9",]$hjust <- 1.5

gg <- ggmap(orlando)
gg <- gg + geom_text(data=median_scores, aes(x=lon, y=lat, label=gsub("School ", "", school)), 
                     hjust=median_scores$hjust, size=3, face="bold", color="darkblue")
gg <- gg + coord_map()
gg <- gg + labs(x=NULL, y=NULL, title=NULL)
gg <- gg + theme_map(base_family="Helvetica")
gg_med_map <- gg

grid.arrange(gg_med_map, gg_med, ncol=2)

Adjust the school labels on the map as needed.

That should help show whatever geographic causality (or lack of) you're trying to show.

Question:

I have a ggplot map, for example:

library(ggmap)
ggmap(get_map())

I'd like the axis labels to be automatically labeled as N-S / W-E: in the above case, for example, instead of lon -95.4 it should show 95.4°E.

I have tried to mess with the scales package and using scale_x_continuous and scale_y_continuous labels and breaks options, but I have not managed to make it work.

It would be awesome to have a scale_y_latitude and scale_x_longitude.

EDIT: Thanks to @Jaap 's answer I got to the following:

scale_x_longitude <- function(xmin=-180, xmax=180, step=1, ...) {
    ewbrks <- seq(xmin,xmax,step)
    ewlbls <- unlist(lapply(ewbrks, function(x) ifelse(x < 0, paste(x, "W"), ifelse(x > 0, paste(x, "E"),x))))
    return(scale_x_continuous("Longitude", breaks = ewbrks, labels = ewlbls, expand = c(0, 0), ...))
}
scale_y_latitude <- function(ymin=-90, ymax=90, step=0.5, ...) {
    nsbrks <- seq(ymin,ymax,step)
    nslbls <- unlist(lapply(nsbrks, function(x) ifelse(x < 0, paste(x, "S"), ifelse(x > 0, paste(x, "N"),x))))
    return(scale_y_continuous("Latitude", breaks = nsbrks, labels = nslbls, expand = c(0, 0), ...))
}

Which works pretty well. But for some reason my R doesn't seem to like the degree symbol in front of the cardinal point... It is displayed as a simple dot, e.g. longitude -24 becomes 24..W


Answer:

Unfortunately, there is no such thing as scale_x_longitude or scale_y_latitude yet. In the meantime here is a workaround in which you specify the labels beforehand:

# load the needed libraries
library(ggplot2)
library(ggmap)

# get the map
m <- get_map(location=c(lon=0,lat=0),zoom=5)

# create the breaks- and label vectors
ewbrks <- seq(-10,10,5)
nsbrks <- seq(-10,10,5)
ewlbls <- unlist(lapply(ewbrks, function(x) ifelse(x < 0, paste(x, "°E"), ifelse(x > 0, paste(x, "°W"),x))))
nslbls <- unlist(lapply(nsbrks, function(x) ifelse(x < 0, paste(x, "°S"), ifelse(x > 0, paste(x, "°N"),x))))

# create the map
ggmap(m) +
  geom_blank() +
  scale_x_continuous(breaks = ewbrks, labels = ewlbls, expand = c(0, 0)) +
  scale_y_continuous(breaks = nsbrks, labels = nslbls, expand = c(0, 0)) +
  theme(axis.text = element_text(size=12))

which gives:


To get the degrees in the functions, you can raise the o as superscript (which will circumvent the need for a special symbol):

scale_x_longitude <- function(xmin=-180, xmax=180, step=1, ...) {
  xbreaks <- seq(xmin,xmax,step)
  xlabels <- unlist(lapply(xbreaks, function(x) ifelse(x < 0, parse(text=paste0(x,"^o", "*W")), ifelse(x > 0, parse(text=paste0(x,"^o", "*E")),x))))
  return(scale_x_continuous("Longitude", breaks = xbreaks, labels = xlabels, expand = c(0, 0), ...))
}
scale_y_latitude <- function(ymin=-90, ymax=90, step=0.5, ...) {
  ybreaks <- seq(ymin,ymax,step)
  ylabels <- unlist(lapply(ybreaks, function(x) ifelse(x < 0, parse(text=paste0(x,"^o", "*S")), ifelse(x > 0, parse(text=paste0(x,"^o", "*N")),x))))
  return(scale_y_continuous("Latitude", breaks = ybreaks, labels = ylabels, expand = c(0, 0), ...))
}    

ggmap(m) +
  geom_blank() +
  scale_x_longitude(xmin=-10, xmax=10, step=5) +
  scale_y_latitude(ymin=-10, ymax=10, step=5) +
  theme(axis.text = element_text(size=12))

which gives the following map:

I used geom_blank just to illustrate the desired effect. You can off course use other geom's (e.g. geom_point) to plot your data on the map.

Question:

I am trying to use GGmap to create a plot of vehicle car crashes by state. The map will have dots which are sized based on the number of car crashes in the state.

In particular I am trying to recreate the usa-plot shown in the visualizing clusters section of this blog post.

However, whenever I try to create the map I get this error.

Error in aperm.default(map, c(2, 1, 3)) : 
  invalid first argument, must be an array

I have setup the Google API and see that it is recieving hits. I have also enabled it and have the key.

In addition I have installed GGmap from the github account using this command:

devtools::install_github("dkahle/ggmap", ref = "tidyup", force=TRUE)

since the CRAN one isn't updated.

I have restarted and quit R several times as well but the error persists.

Even if I just simply run:

get_map()

it still results in the error:

Error in aperm.default(map, c(2, 1, 3)) : 
      invalid first argument, must be an array

Below is my code, it is similar to the code in the blog post:

mydata$State <- as.character(mydata$State)
mydata$MV.Number = as.numeric(mydata$MV.Number)
mydata = mydata[mydata$State != "Alaska", ]
mydata = mydata[mydata$State != "Hawaii", ]
devtools::install_github("dkahle/ggmap", ref = "tidyup", force=TRUE)
library(ggmap)
ggmap::register_google(key = "...") #my key is here
for (i in 1:nrow(mydata)) {
  latlon = geocode(mydata[i,1])
  mydata$lon[i] = as.numeric(latlon[1])
  mydata$lat[i] = as.numeric(latlon[2])
}
mv_num_collisions = data.frame(mydata$MV.Number, mydata$lon, mydata$lat)

colnames(mv_num_collisions) = c('collisions','lon','lat')
usa_center = as.numeric(geocode("United States"))


USAMap = ggmap(get_googlemap(center=usa_center, scale=2, zoom=4), 
extent="normal")
USAMap + 
   geom_point(aes(x=lon, y=lat), data=mv_num_collisions, col="orange", 
alpha=0.4, size=mv_num_collisions$collisions*circle_scale_amt) +  
   scale_size_continuous(range=range(mv_num_collisions$collisions))

I expect the map to output like this

But I cannot seem to get passed this error.

If anyone can help that would be great.

Please let me know if you need any more information.

Thank you.


Answer:

This error is due to the google key not having the appropriate API activity enabled for that key.

Go into the google API console and enable the API "Maps Static API" and it should work for you.

Question:

I can plot the state of Louisiana just fine...

require(ggplot2)
require(ggmap)
require(maps)
LA <- map_data("state", region="louisiana")
ggplot(LA, aes(x=long, y=lat))+geom_polygon()

Now, I have data on how many sales calls were made to particular cities in LA. How would I add a point for each city where a sales call was made to the plot?

salesCalls <- data.frame(State=rep("louisiana",5), 
                             City=c("Baton Rouge", "New Orleans", "Shreveport", "Lafayette", "Mandeville"),
                             Calls=c(10,5,8,13,2))
salesCalls
      State        City Calls
1 louisiana Baton Rouge    10
2 louisiana New Orleans     5
3 louisiana  Shreveport     8
4 louisiana   Lafayette    13
5 louisiana  Mandeville     2

Answer:

require(ggplot2)
require(ggmap)
require(maps)
LA <- map_data("state", region="louisiana")

salesCalls <- data.frame(State=rep("louisiana",5), 
                         City=c("Baton Rouge", "New Orleans", "Shreveport", 
                                "Lafayette", "Mandeville"),
                         Calls=c(10,5,8,13,2))

salesCalls <- cbind(geocode(as.character(salesCalls$City)), salesCalls)

salesCalls
#         lon      lat     State        City Calls
# 1 -91.14032 30.45828 louisiana Baton Rouge    10
# 2 -90.07153 29.95107 louisiana New Orleans     5
# 3 -93.75018 32.52515 louisiana  Shreveport     8
# 4 -92.01984 30.22409 louisiana   Lafayette    13
# 5 -90.06563 30.35825 louisiana  Mandeville     2

ggplot(LA, aes(x=long, y=lat)) +
  geom_polygon() +
  coord_map() +
  geom_point(data=salesCalls, aes(x=lon, y=lat, size=Calls), color="orange")

On a Google Map:

ggmap(get_map(location = 'Louisiana', zoom = 7)) +
  geom_point(data=salesCalls, aes(x=lon, y=lat, size=Calls), color="orange")

Question:

I want to plot a map with a raster overlaying a GoogleMaps base map in ggplot2. Therefore, I used get_map() and insert_raster() like this:

library(ggplot2)
library(ggmap)

bm <- ggmap(get_map(location = "Bangkok", maptype = "hybrid"))

bm + inset_raster(as.raster(r), xmin = r@extent[1], xmax = r@extent[2],
                  ymin = r@extent[3], ymax = r@extent[4])

Is there any possibility to set a alpha and change the fill color?

The result looks like this:


Answer:

Even Faster without fortify:

read the original post below for further information

From this blog entry I found that we can use spatial polygons directly in ggplot::geom_polygon()

r <- raster(system.file("external/test.grd", package="raster"))
# just to make it reproducible with ggmap we have to transform to wgs84
r <- projectRaster(r, crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

rtp <- rasterToPolygons(r)

bm <- ggmap(get_map(location = bbox(rtp), maptype = "hybrid", zoom = 13))
bm + 
  geom_polygon(data = rtp, 
               aes(x = long, y = lat, group = group, 
                   fill = rep(rtp$test, each = 5)), 
               size = 0, 
               alpha = 0.5)  + 
  scale_fill_gradientn("RasterValues", colors = topo.colors(255)) 
How to tackle plotting SPEED if you just need to visualize something

As described below, such plotting might become very slow with large numbers of pixels. Therefore, you might consider to reduce the number of pixels (which in most cases does not really decrease the amount of information in the map) before converting it to polygons. Therefore, raster::aggregate can be used to reduce the number of pixels to a reasonable amount.

The example shows how the number of pixels is decreased by an order of 4 (i.e. 2 * 2, horizontally * vertically). For further information see ?raster::aggregate.

r <- aggregate(r, fact = 2)
#  afterwards continue with rasterToPolygons(r)...
Original Post:

After a while, I found a way to solve this problem. Converting the raster to polygons! This idea then basically was implemented after Marc Needham's blog post.

Yet, there is one drawback: ggplot gets really slow with large numbers of polygons, which you will inevitably face. However, you can speed things up by plotting into a png() (or other) device.


Here is a code example:

library(raster)
library(ggplot2)
library(ggmap)

r <- raster(....) # any raster you want to plot
rtp <- rasterToPolygons(r)
rtp@data$id <- 1:nrow(rtp@data)   # add id column for join

rtpFort <- fortify(rtp, data = rtp@data)
rtpFortMer <- merge(rtpFort, rtp@data, by.x = 'id', by.y = 'id')  # join data

bm <- ggmap(get_map(location = "Shanghai", maptype = "hybrid", zoom = 10))

bm + geom_polygon(data = rtpFortMer, 
                  aes(x = long, y = lat, group = group, fill = layer), 
                  alpha = 0.5, 
                  size = 0) +  ## size = 0 to remove the polygon outlines
     scale_fill_gradientn(colours = topo.colors(255))

This results in something like this:

Question:

I am trying to use ggmap to create a map of the protected areas I am working in with a satellite image from google earth underneath. I can make a very satisfactory image except that it lacks a north arrow and scale bar:

I'm aware that there are very long winded ways to add these elements (e.g. here) but there must surely be a more parsimonious way to do it!

I've tried using map.scale and north.arrow but these both give me:

Error in polygon(xb + arrow.x * s, yb + arrow.y * s, ...) : 
  plot.new has not been called yet

I can get both map.scale and north.arrow to work in base R using plot but then I can't get my satellite image to plot properly. I can also get what I want using arrows and text in base R but again these won't work in ggmap.

The code I'm using is below. You won't have the polygon (so I won't include it in the code) but you'll be able to load the google earth image and replicate the error.

library(rgdal)
library(ggmap)
library(GISTools)

# Load satellite picture

map.centre <- c(lon = 35, lat = -2.5)
map <- get_map(location=map.centre, source="google", maptype="satellite", zoom = 8)

# Plot map

ggmap(map, extent= "device")

map.scale(xc= 34, yc= -3, len= 10, units= "Kilometers",
 ndivs= 4, tcol= "black", scol= "black", sfcol="black")

north.arrow(xb= 35.5, yb= -1, len=100, lab="N")

From doing a bit of reading it seems like the map.scale and north.arrow functions aren't recognising the window the ggmap function creates as an open graphical window. I've done some research and tried to fix this but nothing has worked. Is anyone able to suggest a way to either fix the error I'm receiving or to get a scale bar and north arrow in ggmap without using hundreds of lines of code?


Answer:

It looks like map.scale and north.arrow are designed to work with base graphics, but ggplot uses grid graphics. I'm not that familiar with graphing spatial data, but as a quick hack for the North arrow, the code below includes two different options:

ggmap(map, extent= "device") +
  geom_segment(arrow=arrow(length=unit(3,"mm")), aes(x=33.5,xend=33.5,y=-2.9,yend=-2.6), 
               colour="yellow") +
  annotate(x=33.5, y=-3, label="N", colour="yellow", geom="text", size=4) +
  geom_segment(arrow=arrow(length=unit(4,"mm"), type="closed", angle=40), 
               aes(x=33.7,xend=33.7,y=-2.7,yend=-2.6), colour=hcl(240,50,80)) +
  geom_label(aes(x=33.7, y=-2.75, label="N"),
             size=3, label.padding=unit(1,"mm"), label.r=unit(0.4,"lines"))  

Question:

I'm unsure why none of my data points show up on the map.

   Store_ID visits CRIND_CC  ISCC  EBITDAR top_bottom   Latitude  Longitude
      (int)  (int)    (int) (int)    (dbl)      (chr)     (fctr)     (fctr)
1        92    348    14819 39013 76449.15        top  41.731373  -93.58184
2      2035    289    15584 35961 72454.42        top  41.589428  -93.80785
3        50    266    14117 27262 49775.02        top  41.559017  -93.77287
4       156    266     7797 25095 28645.95        top    41.6143 -93.834404
5        66    234     8314 18718 46325.12        top    41.6002 -93.779236
6       207     18     2159 17999 20097.99     bottom  41.636208 -93.531876
7        59     23    10547 28806 52168.07     bottom   41.56153  -93.88083
8       101     23     1469 11611  7325.45     bottom   41.20982  -93.84298
9       130     26     2670 13561 14348.98     bottom  41.614517  -93.65789
10      130     26     2670 13561 14348.98     bottom 41.6145172  -93.65789
11       24     27    17916 41721 69991.10     bottom  41.597134  -93.49263

> dput(droplevels(top_bottom))
structure(list(Store_ID = c(92L, 2035L, 50L, 156L, 66L, 207L, 
59L, 101L, 130L, 130L, 24L), visits = c(348L, 289L, 266L, 266L, 
234L, 18L, 23L, 23L, 26L, 26L, 27L), CRIND_CC = c(14819L, 15584L, 
14117L, 7797L, 8314L, 2159L, 10547L, 1469L, 2670L, 2670L, 17916L
), ISCC = c(39013L, 35961L, 27262L, 25095L, 18718L, 17999L, 28806L, 
11611L, 13561L, 13561L, 41721L), EBITDAR = c(76449.15, 72454.42, 
49775.02, 28645.95, 46325.12, 20097.99, 52168.07, 7325.45, 14348.98, 
14348.98, 69991.1), top_bottom = c("top", "top", "top", "top", 
"top", "bottom", "bottom", "bottom", "bottom", "bottom", "bottom"
), Latitude = structure(c(11L, 4L, 2L, 7L, 6L, 10L, 3L, 1L, 8L, 
9L, 5L), .Label = c("41.20982", "41.559017", "41.56153", "41.589428", 
"41.597134", "41.6002", "41.6143", "41.614517", "41.6145172", 
"41.636208", "41.731373"), class = "factor"), Longitude = structure(c(3L, 
7L, 5L, 8L, 6L, 2L, 10L, 9L, 4L, 4L, 1L), .Label = c("-93.49263", 
"-93.531876", "-93.58184", "-93.65789", "-93.77287", "-93.779236", 
"-93.80785", "-93.834404", "-93.84298", "-93.88083"), class = "factor")), row.names = c(NA, 
-11L), .Names = c("Store_ID", "visits", "CRIND_CC", "ISCC", "EBITDAR", 
"top_bottom", "Latitude", "Longitude"), class = c("tbl_df", "tbl", 
"data.frame"))

Creating the plot:

map <- qmap('Des Moines') +
       geom_point(data = top_bottom, aes(x = as.numeric(Longitude),
                  y = as.numeric(Latitude)), colour = top_bottom, size = 3)

I get the warning message:

Removed 11 rows containing missing values (geom_point). 

However, this works without the use of ggmap():

ggplot(top_bottom) +  
geom_point(aes(x = as.numeric(Longitude), y = as.numeric(Latitude)),
           colour = top_bottom, size = 3)

How do I get the points to overlay on ggmap??


Answer:

You are using as.numeric() with a factor. As seen here that gives you a level number for the factor (not the number represented). Unsurprisingly, all those levels are points not on the canvas displayed for "Des Moines".

Use as.numeric(as.character(Latitude)) and as.numeric(as.character(Longitude)), as ugly as it seems.

Question:

I am having clipping problems when I try to combine ggmap with shape files. The example in Kahle and Wickham (2013: 158) works fine because the raster image from ggmap covers the entire shape file. Below is an example of what happens when I try to plot the shape file for U.S. states on a ggmap plot that covers a smaller area. The ggmap shows New York City and I want to overlay it with the borders for U.S. states (just as an example). The resulting map doesn't make any sense. The problem is that the shape file gets clipped and ggplot connects the unclipped points. Below is the code. The shape file is from here. I am just showing the last plot here.

How can I solve this problem?

path <- "PATH TO SHAPEFILE"
library("ggmap")
library("rgdal")

# shapefile
states <- readOGR(dsn = path, layer = "states")
states_df <- fortify(states)
# plot shapefile
plot(states, lwd = 0.1)
ggplot(states_df, aes(long, lat, group = group)) +
    geom_polygon(colour = "black", fill = NA, size = 0.1)


# combine ggmap with shapefile
map <- get_map("new york city", zoom = 10, source = "stamen")
ggmap(map, extent = "device")

ggmap(map, extent = "device") +
    geom_polygon(aes(long, lat, group=group), data = states_df, colour = "red", fill = NA, size = 1)

Kahle, David and Hadley Wickham. 2013. “Ggmap: Spatial Visualization with ggplot2.” The R Journal 5(1):144–61.


Answer:

Here is my attempt. I often use GADM shapefiles, which you can directly import using the raster package. I subsetted the shape file for NY, NJ and CT. You may not have to do this in the end, but it is probably better to reduce the amount of data. When I drew the map, ggplot automatically removed data points which stay outside of the bbox of the ggmap image. Therefore, I did not have to do any additional work. I am not sure which shapefile you used. But, GADM's data seem to work well with ggmap images. Hope this helps you.

library(raster)
library(rgdal)
library(rgeos)
library(ggplot2)

### Get data (shapefile)
us <- getData("GADM", country = "US", level = 1)

### Select NY and NJ
states <- subset(us, NAME_1 %in% c("New York", "New Jersey", "Connecticut"))

### SPDF to DF
map <- fortify(states)

## Get a map
mymap <- get_map("new york city", zoom = 10, source = "stamen")


ggmap(mymap) +
geom_map(data = map, map = map, aes(x = long, y = lat, map_id = id, group = group))

If you just want lines, the following would be what you are after.

ggmap(mymap) +
geom_path(data = map, aes(x = long, y = lat, group = group))

Question:

I would like to make a map of the Netherlands with curved lines between cities. I have two dataframes number one called: df_verticles which contains 24 cities with their lat/lon combination. The second dataframe called: df i want to use to draw a curved line between lat/lon from combination to the lat/lon to combination.

> head(df_vertices)
        city AmountSessions   totalkWh AmountRFID scaledAmount   scaledkWh Latitude Longitude
1    Alkmaar          13608  104554.68       1326   0.07139012 0.026941910 52.63903  4.755538
2     Almere          11281  100841.42        930   0.05006999 0.025985067 52.39447  5.282043
3 Amersfoort           7719   67663.30       1198   0.06449876 0.017435647 52.15108  5.383069
4 Amstelveen          25794  236437.93       2616   0.14084204 0.060925915 52.31724  4.859266
5  Amsterdam         402365 3880744.86      18574   1.00000000 1.000000000 52.34560  4.808834

> head(df)
CityChargeSessions   NextCity Amount    sumkWh scaledAmount  scaledkWh Latitude_from Longitude_from Latitude_to Longitude_to
1          Amsterdam    Alkmaar   1058  8133.736   0.18438480 0.15480933      52.34560       4.808834    52.63903     4.755538
2          Amsterdam     Almere    998  7254.133   0.17392820 0.13806786      52.34560       4.808834    52.39447     5.282043
3          Amsterdam Amersfoort    566  4977.404   0.09864064 0.09473489      52.34560       4.808834    52.15108     5.383069
4          Amsterdam Amstelveen   3724 24210.289   0.64900662 0.46079423      52.34560       4.808834    52.31724     4.859266
5             Almere  Amsterdam   1034  8685.526   0.18020216 0.16531155      52.39447       5.282043    52.34560     4.808834
6         Amersfoort  Amsterdam    579  4936.823   0.10090624 0.09396251      52.15108       5.383069    52.34560     4.808834

Normally ggmap is just an underlayer of a ggplot so i decided to plot my desired plot first in ggplot:

ggplot_curve<- ggplot()+
geom_text(data=df_vertices, aes(x = Longitude, y = Latitude+0.025, label = df_vertices$city), size=6)+
geom_point(aes(x = Longitude, y = Latitude, size=scaledkWh), colour="red", data = df_vertices, alpha =0.5)+
scale_size_continuous(range=c(1,30))+
geom_curve(data=df, aes(x=Longitude_from, y=Latitude_from, xend=Longitude_to, yend=Latitude_to),
                    arrow=arrow(angle=15,ends="first",length=unit(0.7,"cm"),type="closed"), size= df$scaledAmount,alpha=0.5, curvature = 0.15)+
theme_bw()

This is exactly what i want, now i just want to add the desired ggmap as underlayer.

Now i just replace ggplot() by ggmap

ggmap_with_curve<- ggmap(map)+
geom_point(aes(x = Longitude, y = Latitude, size=scaledkWh), colour="red", data = df_vertices, alpha =0.5)+
scale_size_continuous(range=c(1,30))+
 geom_curve(data=df, aes(x=Longitude_from, y=Latitude_from, xend=Longitude_to, yend=Latitude_to),
         arrow=arrow(angle=15,ends="first",length=unit(0.7,"cm"),type="closed"), size= df$scaledAmount,alpha=0.5, curvature = 0.15)

As you can see this is not the desired output i was hoping for, R gave me this error message:

geom_curve is not implemented for non-linear coordinates.

I tried to google it, but i coudnt fix it myself.

So my question how do i get this ggplot output with the desired map as underlayer. Could anyone help me out?


Answer:

I think using coord_cartesian() give you what you want. Using the 5 rows showed in your post

map <- get_map(location = "Amsterdam", zoom = 11)
ggmap(map) +
  geom_point(data = df_vertices,
             aes(x = Longitude, y = Latitude, size = scaledkWh),
             colour = "red", alpha =0.5) +
 geom_curve(data = df,
            aes(x = Longitude_from, y = Latitude_from, xend = Longitude_to, yend = Latitude_to),    
            arrow = arrow(angle = 15, ends = "first", length = unit(0.5, "cm"), type = "closed"),
            size = df$scaledAmount, alpha = 0.5, curvature = 0.15, inherit.aes = TRUE)
scale_size_continuous(range=c(1,30)) +
coord_cartesian()

Question:

I'm trying to plot a US map where each state is shaded by the count that it has. I've gotten the shading to work just fine. The problem I'm running into, however, is that the polygons look very jagged (I'm assuming something happened when I tried to merge the map_data('state') with my data frame of counts per state). My data frame before merging has 49 rows (Nevada was missing data in my set), and after merging has many more rows (expected for the long/lat items for the states) but the data appears to be copied correctly for each lat/long pair, so I'm unsure why the poly's are so jagged.

Code:

ggplot() +
  geom_polygon(data=try1, aes(x=long, y=lat, group = group, fill= COUNT)) +
  scale_fill_continuous(low='thistle2', high='darkred', guide='colorbar') +
  theme_bw() + labs(fill="State Map Try Title1", title='Title2', x='', y='') +
  scale_y_continuous(breaks=c()) +
  scale_x_continuous(breaks=c()) +
  theme(panel.border = element_blank())

Any help would be greatly appreciated (and obviously, if there is a better way to do it, I'm open to suggestions!).


Answer:

You don't need to do the merge. You can use geom_map and keep the data separate from the shapes. Here's an example using the built-in USArrests data (reformatted with dplyr):

library(ggplot2)
library(dplyr)

us <- map_data("state")

arr <- USArrests %>% 
  add_rownames("region") %>% 
  mutate(region=tolower(region))

gg <- ggplot()
gg <- gg + geom_map(data=us, map=us,
                    aes(x=long, y=lat, map_id=region),
                    fill="#ffffff", color="#ffffff", size=0.15)
gg <- gg + geom_map(data=arr, map=us,
                    aes(fill=Murder, map_id=region),
                    color="#ffffff", size=0.15)
gg <- gg + scale_fill_continuous(low='thistle2', high='darkred', 
                                 guide='colorbar')
gg <- gg + labs(x=NULL, y=NULL)
gg <- gg + coord_map("albers", lat0 = 39, lat1 = 45) 
gg <- gg + theme(panel.border = element_blank())
gg <- gg + theme(panel.background = element_blank())
gg <- gg + theme(axis.ticks = element_blank())
gg <- gg + theme(axis.text = element_blank())
gg

Question:

With ggmap and ggplot and the following code... (non-reproducible, but imho not necessary to understand the problem).

map <- get_googlemap(center = c(lon = 10.64, lat = 50.56), maptype = "terrain", source = "google", zoom = 6, language = "de-DE", color = "bw")

ggmap(map) + 
  geom_point(data = frage_3_daten, aes(x = lng_google, y = lat_google, colour = pronunciation_id), alpha = 0.2) + 
  scale_colour_hue(name = "Aussprache", labels = c("Krampus", "Grittibänz")) +
  ggtitle("Gebäck in Form einer menschlichen Gestalt") +
  xlab("Länge") + ylab("Breite") +
  theme_srf()

I can produce this beautiful point map on top of the German-speaking Europe.

Now: My only (and hopefully simple) question is: How can I lower the opacity of the background layer, so that the points become more important?

I managed the following "hack" by setting the darken parameter: ggmap(map, darken = c(0.6, "white")).

This almost solves my problem, but maybe there's actually a (hidden) option to globally lower the opacity of the first, map layer (or more generally, of any layer in a plot).


Answer:

If you don't want to adjust the darken parameter, and you don't want to do a deep dive into custom styles for the Google Map, you can modify the ggmap object directly.

The ggmap is essentially a character matrix, where each cell is a hex code for the color to be reproduced there. (There are some extra attributes that describe the longitude and latitude (in the EPSG:4326 coordinate reference system) of the lower left and upper right points of the ggmap, as well as the source and zoom level.)

You can use the adjustcolor() function from base R to take a color (as a hex code, an integer, or a character string-- anything that the col2rgb() function will accept) and dial up or down the red, green, blue, and alpha channels. The alpha channel controls the transparency, where 1 is fully opaque and 0 is fully transparent.

Here's a reproducible example...

First, get the meuse data, which comes with the sp package.

data(meuse)

Transform the dataframe into an sp object, assign it's proper coordinate reference system (which I found here), then transform it's coordinate reference system to longitude/latitude.

coordinates(meuse) = ~x+y
proj4string(meuse) <- "+init=epsg:28992 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs"
meuse <- spTransform(meuse, CRS("+init=epsg:4326"))

Get the Google base map using the approximate center of all the points. I use the same other arguments that you do, but set the zoom much higher since these points are all quite close to each other.

meuse_basemap <- get_map(location = colMeans(coordinates(meuse)), 
                         maptype = "terrain", 
                         source = "google", 
                         zoom = 13, 
                         language = "de-DE", 
                         color = "bw")

Here's where I modify the ggmap to make it transparent. I save the attributes as an object so that I can reassign the same attributes to the more transparent ggmap. I couldn't find a way to modify just the values of the ggmap in place without stripping the ggmap attributes (which then means the modified ggmap doesn't work with the ggmap() function).

meuse_basemap_attributes <- attributes(meuse_basemap)

Create a matrix the same dimensions as the Google basemap, but with all the color hex codes in each cell adjusted to half the transparency.

meuse_basemap_transparent <- matrix(adjustcolor(meuse_basemap, 
                                                alpha.f = 0.5), 
                                    nrow = nrow(meuse_basemap))

Assign the saved attributes to the modified matrix to turn it back into a usable ggmap.

attributes(meuse_basemap_transparent) <- meuse_basemap_attributes

Here's the original plot:

ggmap(meuse_basemap) +
  geom_point(data = as.data.frame(meuse), 
             aes(x = x, y = y, color = dist), 
             cex = 2)

And here's the more transparent one!

ggmap(meuse_basemap_transparent) +
  geom_point(data = as.data.frame(meuse), 
             aes(x = x, y = y, color = dist), 
             cex = 2)

Question:

I'm working with the ggmap package in R and I am relatively new to geospatial data visualizations. I have a data frame of eleven latitude and longitude pairs that I would like to plot on a map, each with a label. Here is the dummy data:

lat<- c(47.597157,47.656322,47.685928,47.752365,47.689297,47.628128,47.627071,47.586349,47.512684,47.571232,47.562283)

lon<-c(-122.312187,-122.318039,-122.31472,-122.345345,-122.377045,-122.370117,-122.368462,-122.331734,-122.294395,-122.33606,-122.379745)

labels<-c("Site 1A","Site 1B","Site 1C","Site 2A","Site 3A","Site 1D","Site 2C","Site 1E","Site 2B","Site 1G","Site 2G")

df<-data.frame(lat,lon,labels)

Now I use annotate to create the data point labels and plot these on a map;

map.data <- get_map(location = c(lon=-122.3485,lat=47.6200), 
                    maptype = 'roadmap', zoom = 11)

pointLabels<-annotate("text",x=uniqueReach$lon,y=c(uniqueReach$lat),size=5,font=3,fontface="bold",family="Helvetica",label=as.vector(uniqueReach$label))

dataPlot <- ggmap(map.data) +
 geom_point(data = uniqueReach,aes(x = df$lon, y = df$lat), alpha = 1,fill="red",pch=21,size = 6) + labs(x = 'Longitude', y = 'Latitude')+pointLabels

This produces a plot of the data points

As you can see, there are two data points that overlap around (-122.44,47.63), and their labels also overlap. Now I can manually add a shift to each label point to keep the labels from overlapping (see this post), but this is not a great technique when I need to produce many of these plots for different sets of latitude and longitude pairs.

Is there a way I can automatically keep data labels from overlapping? I realize whether the labels overlap is dependent on the actual figure size, so I'm open to fixing the figure size at certain dimensions if need be. Thank you in advance for any insights!

EDIT

The following is modified code using the answer given by Sandy Mupratt

# Defining function to draw text boxes
draw.rects.modified <- function(d,...){
  if(is.null(d$box.color))d$box.color <- NA
  if(is.null(d$fill))d$fill <- "grey95"
  for(i in 1:nrow(d)){
    with(d[i,],{
      grid.rect(gp = gpar(col = box.color, fill = fill,alpha=0.7),
                vp = viewport(x, y, w, h, "cm", c(hjust, vjust=0.25), angle=rot))
    })
  }
  d
}


# Defining function to determine text box borders
enlarge.box.modified <- function(d,...){
  if(!"h"%in%names(d))stop("need to have already calculated height and width.")
  calc.borders(within(d,{
    w <- 0.9*w
    h <- 1.1*h
  }))
}

Generating the plot:

dataplot<-ggmap(map.data) + 
                 geom_point(data = df,aes(x = df$lon, y = df$lat), 
                            alpha = 1, fill = "red", pch = 21, size = 6) + 
                  labs(x = 'Longitude', y = 'Latitude') +
                  geom_dl(data = df, 
                      aes(label = labels), 
                      list(dl.trans(y = y + 0.3), "boxes", cex = .8, fontface = "bold"))

This is a MUCH more readable plot, but with one outstanding issue. You'll note that the label "Site 1E" begins to overlap the data point associated with "Site 1A". Does directlabels have a way with dealing with labels overlapping data points belonging to another label?

A final question I have regarding this is how can I plot several duplicate labels using this method. Suppose the labels for data.frame are all the same:

df$labels<-rep("test",dim(df)[1])

When I use the same code, directlabels removes the duplicate label names:

But I want each data point to have a label of "test". Any suggestions?


Answer:

Edit 11 Jan 2016: using ggrepel package with ggplot2 v2.0.0 and ggmap v2.6

ggrepel works well. In the code below, geom_label_repel() shows some of the available parameters.

lat <- c(47.597157,47.656322,47.685928,47.752365,47.689297,47.628128,47.627071,
         47.586349,47.512684,47.571232,47.562283)
lon <- c(-122.312187,-122.318039,-122.31472,-122.345345,-122.377045,-122.370117,
        -122.368462,-122.331734,-122.294395,-122.33606,-122.379745)
labels <- c("Site 1A","Site 1B","Site 1C","Site 2A","Site 3A","Site 1D",
        "Site 2C","Site 1E","Site 2B","Site 1G","Site 2G")

df <- data.frame(lat,lon,labels)

library(ggmap)
library(ggrepel)
library(grid)

map.data <- get_map(location = c(lon = -122.3485, lat = 47.6200), 
                    maptype = 'roadmap', zoom = 11)

ggmap(map.data) + 
   geom_point(data = df, aes(x = lon, y = lat), 
      alpha = 1, fill = "red", pch = 21, size = 5) + 
   labs(x = 'Longitude', y = 'Latitude') +
   geom_label_repel(data = df, aes(x = lon, y = lat, label = labels), 
                 fill = "white", box.padding = unit(.4, "lines"),
                 label.padding = unit(.15, "lines"),
                 segment.color = "red", segment.size = 1)

Original answer but updated for ggplot v2.0.0 and ggmap v2.6

If there is only a small number of overlapping points, then using the "top.bumpup" or "top.bumptwice" method from the direct labels package can separate them. In the code below, I use the geom_dl() function to create and position the labels.

 lat <- c(47.597157,47.656322,47.685928,47.752365,47.689297,47.628128,47.627071,
         47.586349,47.512684,47.571232,47.562283)
 lon <- c(-122.312187,-122.318039,-122.31472,-122.345345,-122.377045,-122.370117,
        -122.368462,-122.331734,-122.294395,-122.33606,-122.379745)
 labels <- c("Site 1A","Site 1B","Site 1C","Site 2A","Site 3A","Site 1D",
        "Site 2C","Site 1E","Site 2B","Site 1G","Site 2G")
 df <- data.frame(lat,lon,labels)

library(ggmap)
library(directlabels)

map.data <- get_map(location = c(lon = -122.3485, lat = 47.6200), 
                    maptype = 'roadmap', zoom = 11)
ggmap(map.data) + 
   geom_point(data = df, aes(x = lon, y = lat), 
      alpha = 1, fill = "red", pch = 21, size = 6) + 
   labs(x = 'Longitude', y = 'Latitude') +
   geom_dl(data = df, aes(label = labels), method = list(dl.trans(y = y + 0.2), 
      "top.bumptwice", cex = .8, fontface = "bold", family = "Helvetica"))

Edit: Adjusting for underlying labels

A couple of methods spring to mind, but neither is entirely satisfactory. But I don't think you will find a solution that will apply to all situations.

Adding a background colour to each label This is a bit of a workaround, but directlabels has a "box" function (i.e., the labels are placed inside a box). It looks like one should be able to modify background fill and border colour in the list in geom_dl, but I can't get it to work. Instead, I take two functions (draw.rects and enlarge.box) from the directlabels website; modify them; and combine the modified functions with the "top.bumptwice" method.

draw.rects.modified <- function(d,...){
  if(is.null(d$box.color))d$box.color <- NA
  if(is.null(d$fill))d$fill <- "grey95"
  for(i in 1:nrow(d)){
    with(d[i,],{
      grid.rect(gp = gpar(col = box.color, fill = fill),
                vp = viewport(x, y, w, h, "cm", c(hjust, vjust=0.25), angle=rot))
    })
  }
  d
}

enlarge.box.modified <- function(d,...){
  if(!"h"%in%names(d))stop("need to have already calculated height and width.")
  calc.borders(within(d,{
    w <- 0.9*w
    h <- 1.1*h
  }))
}

boxes <-
  list("top.bumptwice", "calc.boxes",  "enlarge.box.modified", "draw.rects.modified")

ggmap(map.data) + 
   geom_point(data = df,aes(x = lon, y = lat), 
      alpha = 1, fill = "red", pch = 21, size = 6) + 
   labs(x = 'Longitude', y = 'Latitude') +
   geom_dl(data = df, aes(label = labels), method = list(dl.trans(y = y + 0.3), 
      "boxes", cex = .8, fontface = "bold"))

Add an outline to each label Another option is to use this method to give each label an outline, although it is not immediately clear how it would work with directlabels. Therefore, it would need a manual adjustment of the coordinates, or a search of the dataframe for coordinates that are within a given threshold then adjust. However, here, I use the pointLabel function from maptools package to position the labels. No guarantee that it will work every time, but I got a reasonable result with your data. There is a random element built into it, so you can run it a few time until you get a reasonable result. Also, note that it positions labels in a base plot. The label locations then have to extracted and loaded into the ggplot/ggmap.

lat<- c(47.597157,47.656322,47.685928,47.752365,47.689297,47.628128,47.627071,47.586349,47.512684,47.571232,47.562283)
lon<-c(-122.312187,-122.318039,-122.31472,-122.345345,-122.377045,-122.370117,-122.368462,-122.331734,-122.294395,-122.33606,-122.379745)
labels<-c("Site 1A","Site 1B","Site 1C","Site 2A","Site 3A","Site 1D","Site 2C","Site 1E","Site 2B","Site 1G","Site 2G")
df<-data.frame(lat,lon,labels)

library(ggmap)
library(maptools)  # pointLabel function

# Get map
map.data <- get_map(location = c(lon=-122.3485,lat=47.6200), 
                    maptype = 'roadmap', zoom = 11)

bb = t(attr(map.data, "bb"))   # the map's bounding box

# Base plot to plot points and using pointLabels() to position labels
plot(df$lon, df$lat, pch = 20, cex = 5, col = "red", xlim = bb[c(2,4)], ylim = bb[c(1,3)])
new = pointLabel(df$lon, df$lat, df$labels, pos = 4, offset = 0.5, cex = 1)
new = as.data.frame(new)
new$labels = df$labels

## Draw the map
map = ggmap(map.data) + 
       geom_point(data = df, aes(x = lon, y = lat), 
          alpha = 1, fill = "red", pch = 21, size = 5) + 
       labs(x = 'Longitude', y = 'Latitude') 

## Draw the label outlines 
theta <- seq(pi/16, 2*pi, length.out=32)
xo <- diff(bb[c(2,4)])/400
yo <- diff(bb[c(1,3)])/400

for(i in theta) {
    map <- map + geom_text(data = new,  
       aes_(x = new$x + .01 + cos(i) * xo, y = new$y + sin(i) * yo, label = labels), 
                  size = 3, colour = 'black', vjust = .5, hjust = .8)
}

# Draw the labels
map + 
   geom_text(data = new, aes(x = x + .01, y = y, label=labels), 
     size = 3, colour = 'white', vjust = .5, hjust = .8)

Question:

The goal is to build something like http://rentheatmap.com/sanfrancisco.html

I got map with ggmap and able to plot points on top of it.

library('ggmap')
map <- get_map(location=c(lon=20.46667, lat=44.81667), zoom=12, maptype='roadmap', color='bw')
positions <- data.frame(lon=rnorm(100, mean=20.46667, sd=0.05), lat=rnorm(100, mean=44.81667, sd=0.05), price=rnorm(10, mean=1000, sd=300))
ggmap(map) + geom_point(data=positions, mapping=aes(lon, lat)) + stat_density2d(data=positions, mapping=aes(x=lon, y=lat, fill=..level..), geom="polygon", alpha=0.3)

This is a nice image based on density. Does anybody know how to make something that looks the same, but uses position$property to build contours and scale?

I looked thoroughly through stackoverflow.com and did not find a solution.

EDIT 1

positions$price_cuts <- cut(positions$price, breaks=5)
ggmap(map) + stat_density2d(data=positions, mapping=aes(x=lon, y=lat, fill=price_cuts), alpha=0.3, geom="polygon")

Results in five independent stat_density plots:

EDIT 2 (from hrbrmstr)

positions <- data.frame(lon=rnorm(10000, mean=20.46667, sd=0.05), lat=rnorm(10000, mean=44.81667, sd=0.05), price=rnorm(10, mean=1000, sd=300))
positions$price <- ((20.46667 - positions$lon) ^ 2 + (44.81667 - positions$lat) ^ 2) ^ 0.5 * 10000
positions <- data.frame(lon=rnorm(10000, mean=20.46667, sd=0.05), lat=rnorm(10000, mean=44.81667, sd=0.05))
positions$price <- ((20.46667 - positions$lon) ^ 2 + (44.81667 - positions$lat) ^ 2) ^ 0.5 * 10000
positions <- subset(positions, price < 1000)
positions$price_cuts <- cut(positions$price, breaks=5)
ggmap(map) + geom_hex(data=positions, aes(fill=price_cuts), alpha=0.3)

Results in:

It creates a decent picture on real data as well. This is the best result so far. More suggestions are welcome.

EDIT 3: Here is test data and results of a method above:

https://raw.githubusercontent.com/artem-fedosov/share/master/kernel_smoothing_ggplot.csv

test<-read.csv('test.csv')
ggplot(data=test, aes(lon, lat, fill=price_cuts)) + stat_bin2d(, alpha=0.7) + geom_point() + scale_fill_brewer(palette="Blues")

I believe that there should some method that uses other than density kernel to compute proper polygons. It seems that the feature should be in ggplot out of the box, but I cannot find it.

EDIT 4: I appreciate you time and effort to figure out the proper solution to this seemingly not too complicated question. I voted up both your answers as a good approximations to the goal.

I revealed one problem: the data with circles are too artificial and the approaches do not perform that well on read world data.

Paul's approach gave me the plot:

It seems that it captures patterns of the data that is cool.

jazzurro's approage gave me this plot:

It got the patterns as well. However, both of the plots does not seem to be as beautiful as default stat_density2d plot. I will still wait a couple of days to look if some other solution will come up. If not, I will award the bounty to jazzurro as this will be the result I'll stick to use.

There is an open python + google_maps version of required code. May be someone will find inspiration here: https://github.com/jeffkaufman/apartment_prices


Answer:

Here is my approach. The geom_hex approach is nice. When that came out, I really liked it. I still do. Since you asked something more I tried the following. I think my result is similar to one with stat_density2d. But, I could avoid the issues you had. I basically created a shapefile by myself and drew polygons. I subsetted data by price zone (price_cuts) and drew polygons from the edge to zone center. This approach is in the line of EDIT 1 and 2. I think there is still some distance to reach your ultimate goal if you want to draw a map with a large area. But, I hope this will let you move forward. Finally, I would like to say thank you to a couple of SO users who asked great questions related to polygons. I could not come up with this answer without them.

library(dplyr)
library(data.table)
library(ggmap)
library(sp)
library(rgdal)
library(ggplot2)
library(RColorBrewer)


### Data set by the OP
positions <- data.frame(lon=rnorm(10000, mean=20.46667, sd=0.05), lat=rnorm(10000,    mean=44.81667, sd=0.05))

positions$price <- ((20.46667 - positions$lon) ^ 2 + (44.81667 - positions$lat) ^ 2) ^ 0.5 * 10000

positions <- subset(positions, price < 1000)


### Data arrangement
positions$price_cuts <- cut(positions$price, breaks=5)
positions$price_cuts <- as.character(as.integer(positions$price_cuts))

### Create a copy for now
ana <- positions

### Step 1: Get a map
map <- get_map(location=c(lon=20.46667, lat=44.81667), zoom=11, maptype='roadmap', color='bw')

### Step 2: I need to create SpatialPolygonDataFrame using the original data.
### http://stackoverflow.com/questions/25606512/create-polygon-from-points-and-save-as-shapefile
### For each price zone, create a polygon, SpatialPolygonDataFrame, and convert it
### it data.frame for ggplot.

cats <- list()

for(i in unique(ana$price_cuts)){

foo <- ana %>%
       filter(price_cuts == i) %>%
       select(lon, lat)

    ch <- chull(foo)
    coords <- foo[c(ch, ch[1]), ]

    sp_poly <- SpatialPolygons(list(Polygons(list(Polygon(coords)), ID=1)))

    bob <- fortify(sp_poly)

    bob$area <- i

    cats[[i]] <- bob
}

cathy <- as.data.frame(rbindlist(cats))


### Step 3: Draw a map
### The key thing may be that you subet data for each price_cuts and draw
### polygons from outer side given the following link.
### This link was great. This is exactly what I was thinking.
### http://stackoverflow.com/questions/21748852/choropleth-map-in-ggplot-with-polygons-that-have-holes

ggmap(map) +
    geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)),
                 alpha = .3,
                 data = subset(cathy, area == 5))+
    geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)),
                 alpha = .3,
                 data =subset(cathy, area == 4))+
    geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)),
                 alpha = .3,
                 data = subset(cathy, area == 3))+
    geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)),
                 alpha = .3,
                 data = subset(cathy, area == 2))+
    geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)),
                 alpha= .3,
                 data = subset(cathy, area == 1))+
    geom_point(data = ana, aes(x = lon, y = lat), size = 0.3) +                              
    scale_fill_gradientn(colours = brewer.pal(5,"Spectral")) +
    scale_x_continuous(limits = c(20.35, 20.58), expand = c(0, 0)) +
    scale_y_continuous(limits = c(44.71, 44.93), expand = c(0, 0)) +
    guides(fill = guide_legend(title = "Property price zone"))

Question:

I'm assuming the Raster package has what I need... I'm simply wanting to invert the colors in a Raster image.

The actual scenario is this: I want to invert the raster image returned by a ggmap call:

 library(ggmap)
 ggmap(get_stamenmap(maptype = "toner"))

I want to invert the colors to get a white-on-black version of the Stamen Toner map:


Answer:

This inverts the raster object returned by get_stamenmap()

library("ggmap")
m <- get_stamenmap(maptype = "toner")

# invert colors in raster
invert <- function(x) rgb(t(255-col2rgb(x))/255)    
m_inv <- as.raster(apply(m, 2, invert))

# copy attributes from original object
class(m_inv) <- class(m)
attr(m_inv, "bb") <- attr(m, "bb")

ggmap(m_inv)