Drawn from Neil Saunders’ work at http://rpubs.com/neilfws/229230

… changed for Peru; with a number of bug fixes applied

NOTE: the ggmap 2.6.1 has some bugs; this stuff tries to fix the ones I hit

source(file = "https://raw.githubusercontent.com/dkahle/ggmap/master/R/file_drawer.R")
source(file = "https://raw.githubusercontent.com/dkahle/ggmap/master/R/get_stamenmap.R")
get_stamenmap <- function(bbox = c(left = -95.80204, bottom = 29.38048, right = -94.92313, 
    top = 30.14344), zoom = 10, maptype = c("terrain", "terrain-background", 
    "terrain-labels", "terrain-lines", "toner", "toner-2010", "toner-2011", 
    "toner-background", "toner-hybrid", "toner-labels", "toner-lines", "toner-lite", 
    "watercolor"), crop = TRUE, messaging = FALSE, urlonly = FALSE, color = c("color", 
    "bw"), force = FALSE, where = tempdir(), ...) {
    args <- as.list(match.call(expand.dots = TRUE)[-1])
    argsgiven <- names(args)
    if ("bbox" %in% argsgiven) {
        if (!(is.numeric(bbox) && length(bbox) == 4)) {
            stop("bounding box improperly specified.  see ?get_openstreetmap", 
                call. = F)
        }
    }
    if ("zoom" %in% argsgiven) {
        if (!(is.numeric(zoom) && length(zoom) == 1 && zoom == round(zoom) && 
            zoom >= 0 && zoom <= 18)) {
            stop("scale must be a postive integer 0-18, see ?get_stamenmap.", 
                call. = F)
        }
    }
    if ("messaging" %in% argsgiven) 
        stopifnot(is.logical(messaging))
    if ("urlonly" %in% argsgiven) 
        stopifnot(is.logical(urlonly))
    if ("checkargs" %in% argsgiven) {
        .Deprecated(msg = "checkargs argument deprecated, args are always checked after v2.1.")
    }
    maptype <- match.arg(maptype)
    color <- match.arg(color)
    if (is.null(names(bbox))) 
        names(bbox) <- c("left", "bottom", "right", "top")
    if (maptype %in% c("terrain", "terrain-background", "watercolor")) {
        filetype <- "jpg"
        filetype <- "png"  # error in pre-2.4 function
        # see http://stackoverflow.com/questions/23488022/
        # ggmap-stamen-watercolor-png-error
    } else {
        filetype <- "png"
    }
    fourCorners <- expand.grid(lon = c(bbox["left"], bbox["right"]), lat = c(bbox["bottom"], 
        bbox["top"]))
    fourCorners$zoom <- zoom
    row.names(fourCorners) <- c("lowerleft", "lowerright", "upperleft", "upperright")
    fourCornersTiles <- apply(fourCorners, 1, function(v) LonLat2XY(v[1], v[2], 
        v[3]))
    xsNeeded <- Reduce(":", sort(unique(as.numeric(sapply(fourCornersTiles, 
        function(df) df$X)))))
    numXTiles <- length(xsNeeded)
    ysNeeded <- Reduce(":", sort(unique(as.numeric(sapply(fourCornersTiles, 
        function(df) df$Y)))))
    numYTiles <- length(ysNeeded)
    tilesNeeded <- expand.grid(x = xsNeeded, y = ysNeeded)
    if (nrow(tilesNeeded) > 40) {
        message(paste0(nrow(tilesNeeded), " tiles needed, this may take a while ", 
            "(try a smaller zoom)."))
    }
    base_url <- "http://tile.stamen.com/"
    base_url <- paste(base_url, maptype, "/", zoom, sep = "")
    urls <- paste(base_url, apply(tilesNeeded, 1, paste, collapse = "/"), sep = "/")
    urls <- paste(urls, filetype, sep = ".")
    if (messaging) 
        message(length(urls), " tiles required.")
    if (urlonly) 
        return(urls)
    count <- 0
    nTiles <- nrow(tilesNeeded)
    listOfTiles <- lapply(split(tilesNeeded, 1:nrow(tilesNeeded)), function(v) {
        v <- as.numeric(v)
        get_stamenmap_tile(maptype, zoom, v[1], v[2], force = force, messaging = messaging)
    })
    map <- stitch(listOfTiles)
    if (!crop) {
        attr(map, "source") <- "stamen"
        attr(map, "maptype") <- maptype
        attr(map, "zoom") <- zoom
        return(map)
    }
    if (crop) {
        mbbox <- attr(map, "bb")
        size <- 256 * c(length(xsNeeded), length(ysNeeded))
        slon <- seq(mbbox$ll.lon, mbbox$ur.lon, length.out = size[1])
        slat <- vector("double", length = 256 * length(ysNeeded))
        for (k in seq_along(ysNeeded)) {
            slat[(k - 1) * 256 + 1:256] <- sapply(as.list(0:255), function(y) {
                XY2LonLat(X = xsNeeded[1], Y = ysNeeded[k], zoom, x = 0, y = y)$lat
            })
        }
        slat <- rev(slat)
        keep_x_ndcs <- which(bbox["left"] <= slon & slon <= bbox["right"])
        keep_y_ndcs <- sort(size[2] - which(bbox["bottom"] <= slat & slat <= 
            bbox["top"]))
        croppedmap <- map[keep_y_ndcs, keep_x_ndcs]
    }
    croppedmap <- as.raster(croppedmap)
    class(croppedmap) <- c("ggmap", "raster")
    attr(croppedmap, "bb") <- data.frame(ll.lat = bbox["bottom"], ll.lon = bbox["left"], 
        ur.lat = bbox["top"], ur.lon = bbox["right"])
    attr(croppedmap, "source") <- "stamen"
    attr(croppedmap, "maptype") <- maptype
    attr(croppedmap, "zoom") <- zoom
    croppedmap
}

USGS Data

http://earthquake.usgs.gov/earthquakes/search/

library(ggmap)
library(readr)
library(lubridate)

url <- "http://earthquake.usgs.gov/fdsnws/event/1/query.csv?starttime=2010-1-1%2000%3A00%3A00&endtime=2016-11-25%2023%3A59%3A59&maxlatitude=1.889&minlatitude=-18.896&maxlongitude=-68.071&minlongitude=-83.276&minmagnitude=2.5&eventtype=earthquake&orderby=time"
geonet <- read.csv(url, stringsAsFactors = FALSE)
geonet$Date <- as.Date(geonet$time, format = "%Y-%m-%dT%H:%M:%S", tz = "America/Lima")
geonet$Year <- year(geonet$Date)

Mapping by count

library(jpeg)
library(png)
library(plyr)
Peru.map <- get_stamenmap(bbox = c(left = -83.276, bottom = -25.896, right = -60.071, 
    top = 5), zoom = 4, maptype = "terrain")
p <- ggmap(Peru.map)
p <- p + stat_density_2d(bins = 20, geom = "polygon", size = 2, data = geonet, 
    aes(x = longitude, y = latitude, alpha = ..level.., fill = ..level..))
p <- p + scale_fill_gradient(low = "yellow", high = "red", guide = FALSE) + 
    scale_alpha(range = c(0.02, 0.8), guide = FALSE) + xlab("") + ylab("")
p

Mapping by magnitude

p <- ggmap(Peru.map)
p <- p + geom_point(data = subset(geonet, mag >= 4), aes(x = longitude, y = latitude, 
    color = mag), alpha = 0.5) + scale_color_continuous(low = "yellow", high = "red")
p

### Frequency

barplot(table(geonet$Year), main = "Frequency of +2.5 Quakes in Peru since 2010")

Magnitude

plot(geonet$Year, geonet$mag, pch = 21, xlab = "Year", ylab = "Magnitude", main = "Magnitude of +2.5 Peruvian Earthquakes")