R 的 ggmap 套件：繪製地圖與資料分佈圖，空間資料視覺化

使用 Google 地圖的標記（marker）與路徑（path）：

d <- function(x=-95.36, y=29.76, n,r,a){
  round(data.frame(
    lon = jitter(rep(x,n), amount = a),
    lat = jitter(rep(y,n), amount = a)
  ), digits = r)
}
df <- d(n = 50,r = 3,a = .3)
map <- get_googlemap(markers = df, path = df,, scale = 2)
ggmap(map)

以下我們介紹一些進階的用法，首先產生一些測試用的資料：

mu <- c(-95.3632715, 29.7632836)
nDataSets <- sample(4:10,1)
chkpts <- NULL
for(k in 1:nDataSets){
  a <- rnorm(2); b <- rnorm(2);
  si <- 1/3000 * (outer(a,a) + outer(b,b))
  chkpts <- rbind(chkpts,
    cbind(MASS::mvrnorm(rpois(1,50), jitter(mu, .01), si), k))
}
chkpts <- data.frame(chkpts)
names(chkpts) <- c("lon", "lat","class")
chkpts$class <- factor(chkpts$class)
qplot(lon, lat, data = chkpts, colour = class)

用等高線圖畫在地圖上：

ggmap(get_map(maptype = "satellite"), extent = "device") +
stat_density2d(aes(x = lon, y = lat, colour = class), data = chkpts, bins = 5)

將 crime 資料整理一下：

# only violent crimes
violent_crimes <- subset(crime,
    offense != "auto theft" &
    offense != "theft" &
    offense != "burglary"
    )

# rank violent crimes
violent_crimes$offense <-
factor(violent_crimes$offense,
    levels = c("robbery", "aggravated assault",
      "rape", "murder")
    )

# restrict to downtown
violent_crimes <- subset(violent_crimes,
    -95.39681 <= lon & lon <= -95.34188 &
    29.73631 <= lat & lat <=  29.78400
    )

畫等高線圖：

library(grid)
theme_set(theme_bw(16))
HoustonMap <- qmap("houston", zoom = 14, color = "bw")

# a contour plot
HoustonMap +
stat_density2d(aes(x = lon, y = lat, colour = offense),
    size = 3, bins = 2, alpha = 3/4, data = violent_crimes) +
scale_colour_discrete("Offense", labels = c("Robery","Aggravated Assault","Rape","Murder")) +
theme(
    legend.text = element_text(size = 15, vjust = .5),
    legend.title = element_text(size = 15,face="bold"),
    legend.key.size = unit(1.8,"lines")
    )

二維的 histogram：

# 二維的 histogram
HoustonMap +
stat_bin2d(aes(x = lon, y = lat, colour = offense, fill = offense),
    size = .5, bins = 30, alpha = 2/4, data = violent_crimes) +
scale_colour_discrete("Offense",
    labels = c("Robery","Aggravated Assault","Rape","Murder"),
    guide = FALSE) +
scale_fill_discrete("Offense", labels = c("Robery","Aggravated Assault","Rape","Murder")) +
theme(
    legend.text = element_text(size = 15, vjust = .5),
    legend.title = element_text(size = 15,face="bold"),
    legend.key.size = unit(1.8,"lines")
    )

另一種等高線圖：

HoustonMap +
stat_density2d(aes(x = lon, y = lat, fill = ..level.., alpha = ..level..),
    size = 2, bins = 4, data = violent_crimes, geom = "polygon") +
scale_fill_gradient("ViolentnCrimenDensity") +
  scale_alpha(range = c(.4, .75), guide = FALSE) +
guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10))

加上另外一個圖層：

houston <- get_map("houston", zoom = 14)
overlay <- stat_density2d(aes(x = lon, y = lat, fill = ..level.., alpha = ..level..), bins = 4, geom = "polygon", data = violent_crimes)

HoustonMap +
stat_density2d(aes(x = lon, y = lat, fill = ..level.., alpha = ..level..),
    bins = 4, geom = "polygon", data = violent_crimes) +
scale_fill_gradient("ViolentnCrimenDensity") +
scale_alpha(range = c(.4, .75), guide = FALSE) +
guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10)) +
inset(
    grob = ggplotGrob(ggplot() + overlay +
      scale_fill_gradient("ViolentnCrimenDensity") +
      scale_alpha(range = c(.4, .75), guide = FALSE) +
      theme_inset()
      ),
    xmin = attr(houston,"bb")$ll.lon +
    (7/10) * (attr(houston,"bb")$ur.lon - attr(houston,"bb")$ll.lon),
    xmax = Inf,
    ymin = -Inf,
    ymax = attr(houston,"bb")$ll.lat +
    (3/10) * (attr(houston,"bb")$ur.lat - attr(houston,"bb")$ll.lat)
    )

多張等高線圖：

df <- data.frame(
    x = rnorm(10*100, -95.36258, .05),
    y = rnorm(10*100,  29.76196, .05),
    year = rep(paste("year",format(1:10)), each = 100)
    )
for(k in 0:9){
  df$x[1:100 + 100*k] <- df$x[1:100 + 100*k] + sqrt(.05)*cos(2*pi*k/10)
    df$y[1:100 + 100*k] <- df$y[1:100 + 100*k] + sqrt(.05)*sin(2*pi*k/10)
}

ggmap(get_map(),
    base_layer = ggplot(aes(x = x, y = y), data = df)) +
stat_density2d(aes(fill = ..level.., alpha = ..level..),
    bins = 4, geom = "polygon") +
scale_fill_gradient2(low = "white", mid = "orange", high = "red", midpoint = 10) +
  scale_alpha(range = c(.2, .75), guide = FALSE) +
facet_wrap(~ year)

參考資料：ggmap: Spatial Visualization with ggplot2

繼續閱讀： 123