Visualizing the Weather in Japan: A Sweltering Summer Story!

library(jpndistrict) # install GitHub version to get English names in jpnprefs
library(sf)          # spatial polygons
library(ggplot2)     # plotting essentials
library(dplyr)       # data wrangling
library(tidyr)       # more data wrangling
library(purrr)       # even more data wrangling
library(lubridate)   # date wrangling
library(rvest)       # webscraping
library(polite)      # responsible webscraping
library(stringr)     # string wrangling
library(riem)        # airport weather station data
library(jmastats)    # Japan Meteorological Agency data
library(ggrepel)     # labels
library(gganimate)   # animations
library(geofacet)    # geographical facetting
library(gghighlight) # highlight stuff

# font packages
library(extrafont)
loadfonts(device = "win")

tokyo_his_temp <- jsonlite::read_json("../data/temperature.json", simplifyVector = TRUE)

tokyo_weather_df <- tokyo_his_temp %>% 
  set_names(nm = 1876:2018) %>% 
  map(~as.data.frame(.) %>% 
        modify_if(., is.factor, as.character) %>% 
        modify_if(., is.character, as.numeric)) %>% 
  map2_df(., names(.), ~ mutate(., year = .y)) %>% 
  rename(avg_temp = ".") %>% 
  mutate(year = as.numeric(year))

tokyo_weather_df %>% glimpse()

tokyo_weather_df <- tokyo_weather_df %>% 
  filter(year != 2018) %>% 
  group_by(year) %>% 
  mutate(
    date = seq.Date(from = as.Date("1876-06-01"), 
                    by = "day",
                    length = 122),
    date = format(date, "%m/%d")
  ) %>% 
  ungroup()

glimpse(tokyo_weather_df)

# colorbrewer2.org: diverging 8-class palette
cols <- rev(c('#d53e4f','#f46d43','#fdae61','#fee08b','#e6f598','#abdda4','#66c2a5','#3288bd'))

labels <- c("10", "12", "14", "16", "18", "20", "22", "24", "26", "28", "30", "32")

breaks <- c(seq(10, 32, by = 2))

tokyo_weather_df %>% 
  ggplot(aes(x = date, y = year, fill = avg_temp)) +
  geom_tile() +
  scale_fill_gradientn(
    colours = cols,
    labels = labels,
    breaks = breaks,
    limits = c(11.1, 33.2)) +
  guides(fill = guide_colorbar(title = expression("Temperature " ( degree~C)),
                               reverse = FALSE,
                               title.position = "left",
                               label.position = "bottom",
                               nrow = 1)) +
  scale_y_reverse(limits = c(2017, 1876), expand = c(0, 0),
                  breaks = c(1876, seq(1880, 2015, by = 10), 2017)) +
  scale_x_discrete(breaks = c("06/01", "07/01", "08/01", "09/01", "09/30"),
                   labels = c("June 1st", "July 1st", "Aug. 1st", 
                              "Sept. 1st", "Sept. 30th")) +
  labs(title = "Summers in Tokyo are Getting Longer and Hotter (1876-2017)",
       subtitle = glue::glue("
          One Row = One Year, From June 1st to September 30th
          Average Temperature (Celsius)
          "),
       caption = "Data from Toyo Keizai News via Japan Meteorological Agency") +
  theme_minimal() +
  theme(text = element_text(family = "Roboto Condensed", size = 12),
        axis.title = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(3, "cm"),
        plot.margin=unit(c(1,1,1.5,1.2),"cm"))

url <- "http://www.data.jma.go.jp/obd/stats/etrn/view/annually_s.php?prec_no=44&block_no=47662"

session_jma <- url %>% 
  read_html() %>% 
  html_nodes("#tablefix1") %>% 
  .[[1]] %>% 
  html_table(fill = TRUE, header = FALSE)

tokyo_year_avg_temp <- session_jma %>% 
  select(year = X1, avg_temp = X8, avg_high = X9, avg_low = X10, 
         high_temp = X11, low_temp = X12) %>% 
  slice(-c(1, 2, 3)) %>% 
  mutate(avg_temp = avg_temp %>% str_remove("\\]") %>% as.numeric(),
         avg_high = avg_high %>% str_remove("\\]") %>% as.numeric(),
         avg_low = avg_low %>% str_remove("\\]") %>% as.numeric(),
         high_temp = high_temp %>% str_remove("\\]") %>% as.numeric(),
         low_temp = low_temp %>% str_remove("\\]") %>% as.numeric(),
         year = forcats::as_factor(year))

tokyo_year_avg_temp %>% glimpse()

# 10 Class Divergent Rd-Bu: http://colorbrewer2.org/#type=diverging&scheme=RdBu&n=10
temp_cols <- rev(c('#67001f','#b2182b','#d6604d','#f4a582','#fddbc7',
               '#d1e5f0','#92c5de','#4393c3','#2166ac','#053061'))

tokyo_year_avg_temp %>%
  # filter out 1875 and 2018 due to uncertainties in the measurements
  # http://www.data.jma.go.jp/obd/stats/data/mdrr/man/remark.html
  filter(!year %in% c(1875, 2018)) %>% 
  ggplot(aes(x = year, fill = avg_temp)) +
  geom_bar(position = "fill", width = 1) +
  scale_y_continuous(expand = c(0, 0.01)) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_fill_gradientn(colors = temp_cols, "Average Temperature (Celsius)") +
  labs(title = "Tokyo: Annual Average Temperature (1876-2017)",
       subtitle = "One Stripe = One Year, Left: 1876, Right: 2017") +
  theme_void() +
  theme(text = element_text(family = "Roboto Condensed"),
        legend.position = "bottom",
        legend.title = element_text(family = "Roboto Condensed"))

# web scrape busiest Japan airport codes
session <- bow("https://en.wikipedia.org/wiki/List_of_the_busiest_airports_in_Japan")

japan_airports <- scrape(session) %>% 
  html_nodes("table.wikitable:nth-child(8)") %>% 
  .[[1]] %>% 
  html_table()

japan_airports_clean <- japan_airports %>% 
  janitor::clean_names() %>% 
  mutate(city = cityserved %>% iconv(from = "UTF-8", to = "ASCII//TRANSLIT")) %>% 
  select(airport, city, iata_icao) %>% 
  separate(iata_icao, c("IATA", "ICAO"), "\\/") %>% 
  mutate(city = case_when(
    ICAO == "RJAH" ~ "Ibaraki",            # fix Ibaraki Airport listed as Tokyo
    ICAO == "RJOK" ~ "Kochi",              # fix Kochi, Kochi to Kochi
    TRUE ~ city
  )) %>% 
  # filter out Narita and Kansai Airports
  # filter out RJNA (Nagoya) as not working for some reason
  # RJGG is also Nagoya so I'll use that instead
  filter(!ICAO %in% c("RJAA", "RJBB", "RJNA"))

# grab weather data from stations
summer_weather_riem_raw <-  
  map_df(japan_airports_clean$ICAO, riem_measures,
         date_start = "2018-06-01",
         date_end = "2018-08-31")

# calculate daily averages 
sum_air <- summer_weather_riem_raw %>% 
  mutate(time = as_date(valid) %>% ymd()) %>% 
  separate(time, into = c("year", "month", "day"), sep = "-") %>% 
  group_by(month, day, station) %>% 
  summarize(avg_temp = mean(tmpf),
            max_temp = max(tmpf),
            min_temp = min(tmpf),
            avg_dewp = mean(dwpf)) %>% 
  ungroup() %>% 
  # join riem data with airport data
  left_join(japan_airports_clean, by = c("station" = "ICAO")) %>% 
  glimpse()

library(weathermetrics)
# convert from Fahrenheit to Celsius!
sum_air %>% 
  mutate_at(vars(contains("temp")), 
            ~convert_temperature(temperature = ., old_metric = "f", new_metric = "c")) %>% 
  select(contains("temp")) %>% 
  glimpse()

library(comf)

sum_air %>% 
  mutate_at(vars(contains("avg")), 
            ~convert_temperature(temperature = ., old_metric = "f", new_metric = "c")) %>% 
  mutate(humidex = calcHumx(ta = avg_temp, rh = relh)) %>% 
  select(month, day, station, airport, humidex) %>% 
  slice(1:5)

tky_stations_raw <- jmastats::stations %>% 
  filter(area == "東京" & address != "千代田区北の丸公園")

tky_stations <- tky_stations_raw %>% 
  mutate(
    centroid = map(geometry, st_centroid),
    coords = map(centroid, st_coordinates),
    coord_x = map_dbl(coords, 1),
    coord_y = map_dbl(coords, 2)) %>% 
  select(-centroid, - coords)

library(emojifont)
load.fontawesome()

sf_pref13 <- jpn_pref(pref_code = 13) %>% 
    st_simplify(dTolerance = 0.001) %>% 
    mutate(city_code = as.numeric(city_code)) %>% 
    filter(city_code != 13421) %>%  # filter out the Ogasawara Islands (waaayy off the coast)
    st_union() 

# fix "holes" from st_union(): https://uribo.hatenablog.com/entry/2018/08/14/000541
sf_pref13_FIX <- sf_pref13 %>% 
  st_cast("POLYGON") %>% 
  map(~ .x[1]) %>% 
  st_multipolygon() %>% 
  st_sfc(crs = 4326) %>% 
  as.data.frame() %>% 
  mutate(jis_code = "13", prefecture = "Tokyo") %>% 
  magrittr::set_names(c("geometry", "jis_code", "prefecture")) %>% 
  st_as_sf()

# Insert fontawesome icon into a 'label' column
tky_stations <- tky_stations %>% 
  filter(coord_y > 35) %>% 
  mutate(label = fontawesome(c("fa-university")))

sf_pref13 %>% 
  ggplot() +
  geom_sf(fill = "white") +
  coord_sf(xlim = c(138.9, 139.95),
           ylim = c(35.47, 35.9)) +
  geom_label(data = tky_stations, 
             aes(x = coord_x, y = coord_y,
                 label = label),
             family = "fontawesome-webfont") +
  geom_label_repel(data = tky_stations,
             aes(x = coord_x, y = coord_y,
                 label = station_name),
             size = 3, nudge_y = -0.025) +
  theme_minimal() +
  theme(axis.title = element_blank(),
        axis.ticks = element_blank())

jpnprefs <- jpndistrict::jpnprefs

lat <- jpnprefs %>% select(capital_latitude) %>% as_vector()
lon <- jpnprefs %>% select(capital_longitude) %>% as_vector()

japan_stations_coords <- map2(
  lon, lat, 
  ~ nearest_station(longitude = .x, latitude = .y, 
                    geometry = NULL))

j_stat_ref <- japan_stations_coords %>% 
  map(~as.data.frame(.)) %>% 
  reduce(rbind) 

j_stat_ref %>% 
  summarize(avg_distance = mean(distance) %>% units::set_units(km))

j_station_df <- j_stat_ref %>% 
  select(-geometry) %>% 
  left_join(jmastats::stations %>% st_set_geometry(NULL), 
            by = c("station_no", "area", "station_name")) %>% 
  # erase duplicate stations from the joining
  distinct(station_name, .keep_all = TRUE) 

month <- c(6, 7, 8) # June, July, August
block_no <- j_station_df %>% select(block_no) %>% as_vector()

# create a dataframe consisting of every combination of month and station
df <- crossing(block_no, month) 

# jma_collect()
j_sum_weather_raw <- map2(.x = df$block_no, .y = df$month, 
        ~ jma_collect(item = "daily", block_no = .x, year = 2018, month = .y) %>% 
          mutate(block = .x))

# bind rows and select only the columns necessary
j_all_weather_df <- j_sum_weather_raw %>% 
  bind_rows() %>% 
  select(1, 17, 5, 6, 7) %>% 
  magrittr::set_colnames(c("date", "block", 
                         "temperature_average", "temperature_max", "temperature_min")) 

# stations metadata
stations_df <- jmastats::stations %>% 
  select(pref_code, block_no, station_name, station_no, area) %>% 
  st_set_geometry(NULL) # no need for the station coordinate data now

# combine to get prefecture codes for each station
j_temp_stations_df <- j_all_weather_df %>% 
  left_join(stations_df, by = c("block" = "block_no")) %>% 
  # some stations have same coordinates for 2 station types...
  distinct(block, date, .keep_all = TRUE) %>% 
  mutate(pref_code = as.numeric(pref_code)) %>% 
  left_join(jpnprefs %>% mutate(jis_code = as.numeric(jis_code)),
            by = c("pref_code" = "jis_code")) %>% 
  select(-contains("capital_l"))

# Japan spatial polygons from jpndistrict package
sf_ja <- 1:47 %>% # Prefectural codes 1-47
    map(~jpndistrict::jpn_pref(pref_code = ., district = FALSE)) %>% 
    reduce(rbind) %>% 
    st_simplify(dTolerance = 0.01) %>% 
  mutate(pref_code = as.numeric(pref_code))

# combine temp with spatial df
j_temp_map_stations_df <- j_temp_stations_df %>% 
  mutate(pref_code = as.numeric(pref_code)) %>% 
  left_join(sf_ja, by = c("pref_code", "prefecture"))

cols <- c("Okinawa" = "#e41a1c", "Saga" = "#377eb8", 
          "Kagoshima" = "#4daf4a", "Kumamoto" = "#984ea3",
          "Gifu" = "#ff7f00")

j_temp_colors <- j_temp_stations_df %>% 
    # remove the prefecture suffixes from the english names
    mutate(prefecture_en = str_replace(prefecture_en, "\\-.*", ""))

j_temp_colors %>% 
  ggplot() +
  geom_line(aes(x = date, y = temperature_average, 
                color = forcats::fct_reorder(prefecture_en, temperature_average, 
                                             .fun = mean, .desc = TRUE)), 
            size = 0.75) +
  scale_color_manual(name = "Top 5 Prefectures", values = cols) +
  gghighlight(mean(temperature_average), max_highlight = 5,
              use_direct_label = FALSE) +
  scale_x_date(labels = c("June 1st", "July 1st", "August 1st", "August 31st"),
               breaks = as.Date(c("2018-06-01", "2018-07-01", "2018-08-01", "2018-08-31")),
               expand = c(0, 0)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 5),
                     expand = c(0, 0)) +
  labs(title = "Top 5 Prefectures with Highest Average Temperature", 
       subtitle = "Averaged over the Summer Months (June 1st - August 31st)",
       caption = "Source: Japan Meteorological Agency",
       y = "Average Temperature (°C)") +
  theme_minimal() +
  theme(axis.title.x = element_blank(),
        text = element_text("Roboto Condensed"))

region_plots <- j_temp_colors %>% 
  filter(date == "2018-06-21") %>% # day of the summer solstice
  group_by(region_en) %>% 
  mutate(region_avg_temp = mean(temperature_average) %>% round(digits = 2)) %>% 
  nest() %>% 
  mutate(plot = map2(.x = region_en, .y = data, 
                     ~ ggplot(data = .y,
                              aes(
                                x = forcats::fct_reorder(prefecture_en, temperature_average, 
                                                         .desc = TRUE),
                                y = temperature_average
                              )) +
                       geom_segment(aes(y = region_avg_temp, yend = temperature_average,
                                        xend = prefecture_en), 
                                    color = "skyblue",
                                    size = 1.5) +
                       geom_point(color = "red", 
                                  size = 2.5,
                                  show.legend = FALSE) +
                       geom_hline(aes(yintercept = region_avg_temp)) +
                       scale_x_discrete(expand = c(0.01, 0.01)) +
                       scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) +
                       labs(title = glue::glue(
                         "{.x} Region"),
                         subtitle = glue::glue(
                           "Average Regional Temperature ({.y$region_avg_temp} °C)"),
                         x = NULL,
                         y = "Temperature (°C)") +
                       theme_minimal() +
                       theme(text = element_text(family = "Roboto Condensed", size = 7),
                             axis.text.x = element_text(size = 8.5),
                             panel.grid.major.x = element_blank(),
                             panel.grid.minor.y = element_blank()))
  )

glimpse(region_plots)

# show ALL 8 plots at once:
# walk(region_plots$plot, print)

kanto <- region_plots$plot[[1]]
kansai <- region_plots$plot[[2]]
hokkaido <- region_plots$plot[[3]]
tohoku <- region_plots$plot[[4]]
chubu <- region_plots$plot[[5]]
chugoku <- region_plots$plot[[6]]
kyushu <- region_plots$plot[[7]]
shikoku <- region_plots$plot[[8]]

library(patchwork)
# group them up, going top = north, bottom = south
(hokkaido + tohoku) / (chubu + kanto) / (kansai + chugoku) / (shikoku + kyushu) +
  plot_annotation(title = "Regional Temperatures on the Summer Solstice (June 21st, 2018)",
                  subtitle = "Horizontal line: Average Regional Temperature",
                  caption = "Source: Japan Meteorological Agency")

j_temp_map_stations_df %>% 
  ggplot() +
  geom_sf(aes(fill = temperature_average)) +
  scale_fill_gradientn(colours = jmastats:::jma_pal(
    palette = "relative", .attribute = FALSE)[6:1],
    labels = c("35", "30", "25", 
               "20",  "15", "10"),
    breaks = c(35, 30, 25, 20, 15, 10),
    limits = c(9, 35.5),
    name = glue::glue(
      "Temperature (°C)")) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(
    title = "Average Daily Temperature (°C) in Japan",
    subtitle = "{frame_time}",
    caption = "Source: Japan Meteorological Agency",
    x = NULL, y = NULL) +
  theme_minimal() +
  theme(
    text = element_text(family = "Roboto Condensed"),
    legend.position = c(0.8, 0.3),
    legend.background = element_rect(color = "black"),
    legend.title = element_text(size = 10),
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    plot.margin = unit(c(1, 0, 0, 0), "mm")) +
  # animate on date
  transition_time(date, range = as.Date(c("2018-06-01", "2018-08-31")))

# Fix issue of grid dataframe not recognizing the facetted variable name:
pref_names <- j_temp_map_stations_df %>% 
  mutate(prefecture_en = str_replace(prefecture_en, "\\-.*", "")) %>% 
  distinct(pref_code, prefecture_en) %>% 
  arrange(pref_code) %>% 
  magrittr::use_series(prefecture_en)

jp_prefs_grid1 <- jp_prefs_grid1 %>% 
  arrange(code_pref_jis) %>% 
  # fix Tochigi, Gunma & Ibaraki positioning:
  mutate(col = as.numeric(case_when( 
  code == "9" ~ "13",
  code == "8" ~ "14",
  code == "10" ~ "12",
  TRUE ~ as.character(col))),
  row = as.numeric(case_when( 
    code == "9" ~ "5",
    code == "8" ~ "6",
    TRUE ~ as.character(row)))) %>% 
  mutate(prefecture_en = pref_names)

# colors by region
cols <- c("Kyushu" = "#e41a1c", "Shikoku" = "#377eb8", "Chugoku" = "#4daf4a", "Chubu" = "#984ea3",
          "Kansai" = "#ff7f00", "Kanto" = "#ffff33", "Tohoku" = "#a65628", "Hokkaido" = "#f781bf")

# plot!
j_temp_colors %>% 
  ggplot(aes(date, temperature_average)) +
  geom_line(aes(color = region_en), size = 0.75, show.legend = FALSE) +
  scale_color_manual(values = cols) +
  scale_y_continuous(breaks = c(10, 20, 30)) +
  gghighlight(use_direct_label = FALSE) +
  labs(title = glue::glue("Average Daily Temperature (°C) in Japan"),
       subtitle = "June ~ August 2018, Colored by Region",
       caption = "Source: Japan Meteorological Agency") +
  theme_minimal() +
  theme(text = element_text(family = "Roboto Condensed"),
        title = element_text(size = 18), 
        axis.title = element_blank(),
        panel.grid.minor.x = element_blank(),
        strip.background = element_rect(color = "black", linetype = "solid"),
        strip.text.x = element_text(size = 10.5),
        axis.text = element_text(size = 8)) +
  facet_geo(~ prefecture_en, 
            grid = "jp_prefs_grid1")