jma_temperature_plots.R

# Requirements ------------------------------------------------------------
devtools::install_git(
  "https://gitlab.com/uribo/jmastats"
)
remotes::install_github("yutannihilation/gghighlight")
remotes::install_github("ropensci/rnaturalearth")
remotes::install_github("ropensci/rnaturalearthhires")
remotes::install_github("thomasp85/gganimate", ref = "c4e9c9ae2e338589257989f953a75edd76579696")
# R Session ---------------------------------------------------------------
library(jmastats)
library(tidyverse)
library(sf)
quartzFonts(ipa = quartzFont(rep("IPAPGothic", 4)))
theme_set(theme_classic(base_size = 12, base_family = "IPAPGothic"))

# 47都道府県・エリアから1地点ずつ、観測地点のデータを選択する
set.seed(125)
df_target_stations <-
  stations %>%
  # 観測装置の種類が「有線ロボット気象計」(いわゆるアメダス) であるものに制限
  filter(station_type %in% c("四")) %>%
  filter(block_no != "0092") %>%
  group_by(pref_code, area) %>%
  sample_n(1) %>%
  ungroup() %>%
  select(pref_code, area, station_name, block_no)

df_targets <-
    purrr::map2_dfr(
      .x = df_target_stations$block_no %>%
        set_names(),
      .y = df_target_stations$station_name,
      ~ jma_collect(item = "hourly",
                    block_no = .x,
                    year = 2018,
                    month = 7,
                    day = 31) %>%
        select(date, time, starts_with("temperature")) %>%
        parse_unit() %>%
        mutate(station = .y),
      .id = "block_no")
  
df_temperature_summary <-
  df_targets %>%
  group_by(block_no, station) %>%
  summarise_at(.vars = vars(temperature), .funs = funs(min, max), na.rm = TRUE) %>%
  ungroup() %>%
  left_join(stations %>%
              st_set_geometry(NULL) %>% 
              select(block_no, area), by = "block_no")

# jma_124stations_temperature_days20180731 ----------------------------
df_plot <-
  df_temperature_summary %>%
  tidyr::gather(type, value, -block_no, -station, -area)

p_out <- 
  df_plot %>% {
  ggplot(data = ., aes(forcats::fct_reorder(station, value), value)) +
    geom_line(size = 1.2, color = "skyblue") +
    geom_point(aes(color = value), size = 2) +
    scale_color_gradientn(colours = jmastats:::jma_pal(palette = "relative", .attribute = FALSE)[5:1],
                          labels = c("35~", "30~35", "25~30", "20~25",
                                     "15~20"),
                          breaks = c(35, 30, 25, 20, 15),
                          limits = c(15, 37)) +
    geom_text(data = dplyr::filter(., type == "max"), 
              aes(label = value), nudge_x = 0.5, nudge_y = -0.6, size = 2.2) +
    geom_text(data = dplyr::filter(., type == "min"), 
              aes(label = value), nudge_x = 0.5, nudge_y = 0.6, size = 2.2) +
    coord_flip() +
    theme_minimal(base_size = 8) +
    guides(color = guide_legend(title = paste("Temperature", "[\u2103]"),
                                reverse = TRUE,
                                title.position = "left",
                                label.position = "bottom",
                                keywidth = 2,
                                nrow = 1)) +
    theme(legend.position = "bottom") +
    scale_x_discrete(expand = c(0,1)) +
    scale_y_continuous(breaks = seq(0, 40, by = 10),
                       minor_breaks = seq(0, 40, by = 5)) +
    # guides(color = FALSE) +
    labs(x = NULL,
         y = NULL,
         title = "2018年7月31日の各気象観測所における気温",
         subtitle = "日最低・最高気温を表示",
         caption = "Source: 気象庁\nhttp://www.data.jma.go.jp/obd/stats/etrn/index.php")
}

ggsave("jma_124stations_temperature_days20180731.png", p_out, width = 5, height = 8, dpi = "retina")


# 2. jma_124stations_temperature_days20180731_regional_ts ----------------------------------------------------------------------
df_plot2 <-
  df_targets %>%
  left_join(stations %>%
              select(block_no, area, pref_code), by = "block_no") %>%
  left_join(jpndistrict::jpnprefs %>%
              select(pref_code = jis_code, region), by = "pref_code") %>% 
  st_sf()

library(gghighlight)
p_out2 <- 
  df_plot2 %>%
  mutate(time = hms::hms(hours = time),
         region = forcats::fct_inorder(region)) %>%
  ggplot(aes(time, temperature, group = station)) +
  geom_line(aes(color = temperature)) +
  scale_color_gradientn(colours = jmastats:::jma_pal(palette = "relative", .attribute = FALSE)[5:1],
                        labels = c("35~", "30~35", "25~30", "20~25",
                                   "15~20"),
                        breaks = c(35, 30, 25, 20, 15),
                        limits = c(15, 37)) +
  scale_x_time(labels = scales::date_format("%H:%M")) +
  gghighlight(use_direct_label = FALSE, unhighlighted_colour = ggplot2::alpha("grey", 0.3)) +
  facet_wrap(~ region, ncol = 1) +
  theme_minimal() +
  guides(color = guide_legend(title = paste("Temperature", "[\u2103]"),
                              reverse = TRUE,
                              title.position = "left",
                              label.position = "bottom",
                              keywidth = 2,
                              nrow = 1)) +
  labs(x = NULL,
       y = NULL,
       title = "2018年7月31日の地域別での各気象観測所における気温の時系列推移",
       subtitle = "1時間毎の観測値を表示",
       caption = "Source: 気象庁\nhttp://www.data.jma.go.jp/obd/stats/etrn/index.php") +
  theme(rect        = element_blank(), 
        panel.grid  = element_blank(),
        axis.line  = element_blank(),
        axis.text.y = element_blank(),
        legend.position = "bottom",
        strip.text = element_text(angle = 0, hjust = 0.05, vjust = 0.2))

ggsave("jma_124stations_temperature_days20180731_regional_ts.png", p_out2, width = 5, height = 8, dpi = "retina")

# 3. jma_124stations_temperature_days20180731_map ----------------------------------------------------------------------
library(rnaturalearth)
library(ggrepel)
ne_jpn <-
  ne_states(country = "Japan", returnclass = "sf") %>%
  tibble::new_tibble(subclass = "sf") %>% 
  st_union()

ne_jpn <-
  st_crop(
    ne_jpn,
    st_polygon(list(
      rbind(
        st_point(c(123, 24)),
        st_point(c(147, 24)),
        st_point(c(147, 45.5)),
        st_point(c(123, 45.5)),
        st_point(c(123 ,24))
      ))) %>%
      st_sfc(crs = 4326)
  ) %>% 
  st_transform(crs = "+proj=laea +lat_0=30 +lon_0=165")

df_plot3 <-
  df_plot2 %>%
  st_transform(crs = "+proj=laea +lat_0=30 +lon_0=165") %>%
  mutate(time = hms::hms(hours = time),
         longitude = sf::st_coordinates(geometry)[, 1],
         latitude = sf::st_coordinates(geometry)[, 2]) %>%
  st_set_geometry(NULL)

df_plot3_static <- 
  df_plot3 %>% 
  group_by(station) %>% 
  filter(time == hms::as.hms("12:00:00")) %>% 
  ungroup()

p_out3 <- 
  df_plot3_static %>% {
    ggplot() +
      geom_sf(data = ne_jpn) +
      geom_linerange(data = .,
                     aes(longitude,
                         ymin = latitude,
                         ymax = latitude + 50000,
                         color = temperature),
                     size = 1.2,
                     linetype = 1) +
      scale_color_gradientn(colours = jmastats:::jma_pal(palette = "relative", .attribute = FALSE)[5:1],
                            labels = c("35~", "30~35", "25~30", "20~25",
                                       "15~20"),
                            breaks = c(35, 30, 25, 20, 15),
                            limits = c(15, 37)) +
      geom_label_repel(data = dplyr::group_by(., region) %>% 
                  arrange(desc(temperature)) %>% 
                  slice(1L) %>% 
                  ungroup(), 
                aes(x = longitude, y = latitude, label = paste0(station, "\n", temperature)), 
                size = 1.8,
                segment.color = "gray50",
                segment.alpha = 0.3,
                label.padding = unit(0.15, "lines"),
                nudge_x = 10000,
                nudge_y = 90000) +
      geom_segment(aes(x = longitude, y = latitude, xend = longitude + 22000, yend = latitude - 50000),
                   data = .,
                   color = "gray20",
                   alpha = 0.1,
                   size = 1.0) +
      labs(x = NULL,
           y = NULL,
           title = "2018年7月31日正午の気温",
           subtitle = "各地域での最高気温の地点をラベルで表示",
           caption = "Source: 気象庁\nhttp://www.data.jma.go.jp/obd/stats/etrn/index.php") +
      guides(color = guide_legend(title = paste("Temperature", "[\u2103]"),
                                  reverse = TRUE,
                                  title.position = "top",
                                  label.position = "left",
                                  keywidth = 0.5,
                                  nrow = 1)) +
      theme_minimal(base_size = 8) +
      theme(rect = element_blank(),
            axis.line = element_blank(),
            axis.text = element_blank(),
            axis.ticks = element_blank(),
            legend.position = c(0.1, 0.9),
            legend.box.background = element_rect(colour = "grey50"),
            legend.justification = c(0, 1))
  }  

ggsave("jma_124stations_temperature_days20180731_map.png", p_out3, width = 5, height = 5, dpi = "retina")


library(gganimate)
p_out3_animation <- 
  df_plot3 %>% 
  mutate(time = lubridate::ymd_hms(paste(date, as.character(time)))) %>% 
  ggplot() +
  geom_sf(data = ne_jpn) +
  geom_linerange(aes(longitude,
                     ymin = latitude,
                     ymax = latitude + 50000,
                     color = temperature),
                 size = 1.2,
                 linetype = 1) +
  scale_color_gradientn(colours = jmastats:::jma_pal(palette = "relative", .attribute = FALSE)[5:1],
                        labels = c("35~", "30~35", "25~30", "20~25",
                                   "15~20"),
                        breaks = c(35, 30, 25, 20, 15),
                        limits = c(15, 37)) +
  geom_segment(aes(x = longitude, y = latitude, xend = longitude + 22000, yend = latitude - 50000),
               color = "gray20",
               alpha = 0.1,
               size = 1.0) +
  labs(x = NULL,
       y = NULL,
       title = "2018年7月31日の気温",
       subtitle = "時刻: {frame_time}",
       caption = "Source: 気象庁\nhttp://www.data.jma.go.jp/obd/stats/etrn/index.php") +
  guides(color = guide_legend(title = paste("Temperature", "[\u2103]"),
                              reverse = TRUE,
                              title.position = "top",
                              label.position = "left",
                              keywidth = 0.5,
                              nrow = 1)) +
  theme_minimal(base_size = 12) +
  theme(rect = element_blank(),
        axis.line = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        legend.position = c(0.1, 0.9),
        legend.box.background = element_rect(colour = "grey50"),
        legend.justification = c(0, 1)) +
  transition_time(time)

magick::image_write(animate(p_out3_animation, nframes = 80),
                    "jma_124stations_temperature_days20180731_map.gif")