過去７日間の地震マップ

この記事を参考に、ごく基本的な地図作成のスクリプト。

データソースはUSGSの地震データで、過去１時間から１日、７日、３０日、また、マグニチュードの下限についても、いくつか選択できる。

ここでは、Past 7 Days > All Earthquakes を用いる。

ダウンロードした csv ファイルから、日付のフォーマットを修正し、続いて太平洋中心の地図（map_data の world2）にプロットするために、通常の経度：-180°〜180° を 0°〜360° に変換する処理などに dplyr の関数を利用しているのがポイント。

library(lubridate)
library(dplyr)
library(ggplot2)
library(ggthemes)

url <- "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_week.csv"
download.file(url = url, destfile = "all_week.csv")
df <- read.csv("all_week.csv", stringsAsFactors = FALSE)

df$time <- ymd_hms(df$time)
df <- arrange(df, -row_number())
quakes <- df %>% filter(type == "earthquake") %>%
   mutate(longitude = ifelse(longitude < 0, longitude + 360, longitude))

quakes$mag <- factor(round(quakes$mag))

map <- map_data('world2') %>% filter(region != "Antarctica")

title <- paste("Earthquakes map from",
                      as.Date(quakes$time[1]), "to", as.Date(quakes$time[nrow(df)]))

p <- ggplot()
p <- p + geom_map(data = map, map = map,
                                 aes(x = long, y = lat, group = group, map_id = region),
                                 fill = "grey60", colour = "grey90", size = 0.1)
p <- p + geom_point(data = quakes,
                                  aes(x = longitude, y = latitude, colour = mag))
p <- p + scale_color_brewer(palette = 8)
p <- p + theme_map() + coord_equal() + ggtitle(title)
p <- p + theme(legend.position = "bottom")
p <- p + guides(colour = guide_legend(nrow = 1))
print(p)