h3forrお試し

h3forr_sample.R

library(h3forr)
library(jpmesh)
library(sf)
library(dplyr)
library(mapview)
sfc_poly_to_h3 <- function(geometry, res = 7) {
  coords <- sf::st_coordinates(sf::st_centroid(st_geometry(geometry)))
  h3forr::geo_to_h3(c(coords[2], coords[1]), res = res)
}

sf_pref01_mesh10km <-
  administration_mesh(code = 1, to_mesh_size = 10)

sf_pref01_mesh10km_h3 <- 
  sf_pref01_mesh10km %>% 
  mutate(h3 = purrr::pmap_chr(
    .,
    ~ sfc_poly_to_h3(..2, res = 5))) %>% 
  pull(h3) %>% 
  unique() %>% 
  h3_to_geo_boundary() %>% 
  geo_boundary_to_sf()

sf_pref01_mesh80km_h3 <- 
  sf_pref01 %>% 
  mutate(h3 = purrr::pmap_chr(
    .,
    ~ sfc_poly_to_h3(..2, res = 3))) %>% 
  pull(h3) %>% 
  unique() %>% 
  h3_to_geo_boundary() %>% 
  geo_boundary_to_sf()

mapview(sf_pref01_mesh80km_h3, col.regions = "#F4A935") +
  mapview(sf_pref01_mesh10km_h3, col.regions = "#4886A5")

jpn_pref01_landuse.R

###########################################
# h3 (h3forr) による六角形ポリゴンでの地図作成
# データ: 国土交通省 国土数値情報 (土地利用3次メッシュデータ 第2.6版 L03-a 平成26年度 http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L03-a.html )
###########################################
# 1/3 Load pkgs and define function -------------------------------------------------------------------
# remotes::install_github("/crazycapivara/h3forr")
library(h3forr)
library(ggplot2)
library(jpmesh)
library(sf)
library(dplyr)
library(mapview)
# remotes::install_github("uribo/zipangu")
library(zipangu)
# remotes::install_github("thomasp85/patchwork")
library(patchwork)
ksj_l03_tidy <- function(path) {
  zipangu:::read_ksj_l03a(path) %>% 
    sf::st_transform(crs = 4326) %>% 
    tidyr::gather("type", "area", -meshcode, -geometry) %>% 
    dplyr::mutate(area = units::set_units(area, m^2) %>% 
                    units::set_units(km^2) %>% 
                    units::drop_units()) %>% 
    dplyr::select(names(.)[!names(.) %in% attr(., "sf_column")])
}
mesh_to_h3 <- function(data, type, res = 7) {
  data %>%
    dplyr::filter(type %in% !! rlang::enquo(type)) %>%
    dplyr::mutate(h3 = purrr::pmap_chr(
      .,
      function(geometry, ...) {
        coords <- sf::st_coordinates(sf::st_centroid(sf::st_geometry(obj = geometry)))
        h3forr::geo_to_h3(c(coords[2], coords[1]), res = res)
      })) %>%
    sf::st_drop_geometry() %>%
    dplyr::group_by(h3, type) %>%
    dplyr::summarise(area = sum(area, na.rm = TRUE)) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(geomety = h3forr::h3_to_geo_boundary(h3) %>%
                    h3forr::geo_boundary_to_sf() %>%
                    sf::st_geometry()) %>%
    sf::st_sf()
}
summarise_mesh10km_typearea_mesh <- function(data, type) {
  type <- rlang::enquo(type)
  data %>%
    dplyr::filter(type %in% !! type) %>%
    mutate(meshcode = stringr::str_sub(meshcode, 1, 6)) %>%
    group_by(meshcode, type) %>%
    # proporion... dataはdf_pref01_mesh
    # summarise(area = sum(area) / n()) %>%
    summarise(area = sum(area, na.rm = TRUE)) %>% 
    ungroup() %>% 
    meshcode_sf(meshcode)
}
summarise_mesh10km_typearea_h3 <- function(data, type) {
  type <- rlang::enquo(type)
  data %>%
    dplyr::filter(type %in% !! type) %>% 
    dplyr::mutate(meshcode = stringr::str_sub(meshcode, 1, 6)) %>%
    dplyr::group_by(meshcode, type) %>%
    dplyr::summarise(area = sum(area, na.rm = TRUE)) %>%
    dplyr::ungroup() %>%
    mesh_to_h3(type = !! type, res = 5) %>% 
    dplyr::group_by(h3, type) %>% 
    dplyr::summarise(area = sum(area, na.rm = TRUE)) %>%
    dplyr::ungroup()
}
plot_type <- function(data, type) {
  type <- rlang::enquo(type)
  d <- 
    df_landuse_cd_h26 %>%
    dplyr::filter(type %in% !! type)
  p <- 
    data %>% 
    ggplot() +
    geom_sf(aes(fill = area),
            color = "#00000060",
            size = 0.2) +
    coord_sf(datum = NA) +
    theme_bw(base_size = 4, 
             base_family = "IPAexGothic") + 
    theme(rect = element_blank(), 
          text = element_text(family = "IPAexGothic"),
          legend.position = "top",
          legend.direction = "horizontal",
          legend.key.height = unit(0.2, "line"),
          legend.key.width = unit(0.5, "line")) +
    guides(fill = guide_colorbar(title = bquote(.(d$type_jp) ~ (km^2)), 
                                 title.position = "top", 
                                 title.hjust = 0.5))
  col <-
    d %>%
    dplyr::pull("col")
  p + scale_fill_gradientn(colours = c("#FFFFFF", col))
}

# 2/3 Setup dataset --------------------------------------------------------------------
landuse_type <- c("ta", "nouchi", "sinrin", "kouch", "tatemono", "douro", 
                  "tetudou", "other", "kasen", "kaihin", "kaisui", "golf", NA_character_) %>% 
  purrr::set_names(c("田", "その他の農用地", "森林", "荒地", "建物用地", "道路", "鉄道", 
                     "その他の用地", "河川地及び湖沼", "海浜", "海水域", "ゴルフ場", "解析範囲外"))
df_landuse_cd_h26 <- 
  data.frame(code = c(10, 20, 50, 60, 70, 91, 92, 100, 110, 140, 150, 160, 0),
             type = landuse_type,
             col = apply(matrix(
               c(255, 255, 0,
                 255, 204, 153,
                 0, 170, 0,
                 255, 153, 0,
                 255, 0, 0,
                 140, 140, 140,
                 180, 180, 180,
                 200, 70, 15,
                 0, 0, 255,
                 255, 255, 153,
                 0, 204, 255,
                 0, 255, 0,
                 255, 255, 255),
               ncol = 3,
               byrow = TRUE), 1, function(x) {
                 rgb(red = x[1], 
                     green = x[2], 
                     blue = x[3], 
                     maxColorValue = 255)
               }),
             stringsAsFactors = FALSE) %>% 
  tibble::rownames_to_column(var = "type_jp")
landuse_type <- 
  landuse_type %>% 
  purrr::keep(~ !is.na(.x))
sf_pref01_mesh10km <-
  administration_mesh(code = 1, to_mesh_size = 10)
sf_pref01 <- 
  rnaturalearth::ne_states(country = "Japan", returnclass = "sf") %>% 
  tibble::new_tibble(nrow = nrow(.), subclass = "sf") %>% 
  filter(iso_3166_2 == "JP-01") %>% 
  st_geometry() %>% 
  st_sf()

# データをダウンロードする場合
# df_pref01_ksj_l03 <- 
#   c(sf_pref01_mesh10km %>% 
#     pull(meshcode) %>% 
#     stringr::str_sub(1, 4) %>% 
#     unique()) %>% 
#   purrr::map(
#     ~  zipangu::read_ksj_l03a(.year = 2014, .meshcode = .x, .download = TRUE)
#   ) %>% 
#   purrr::reduce(rbind)
df_pref01_ksj_l03 <- 
  fs::dir_ls("~/Documents/resources/国土数値情報/L03-a/", 
             recurse = TRUE, 
             regexp = "-jgd_GML/.+.shp") %>% 
  purrr::map(ksj_l03_tidy) %>% 
  purrr::reduce(rbind)

# 県単位で1kmメッシュは細かいので10kmにして集計、可視化
# summarise_mesh10km_typearea_mesh(df_pref01_ksj_l03, "sinrin")
# summarise_mesh10km_typearea_h3(df_pref01_ksj_l03, "sinrin")
df_plot_mesh <- 
  landuse_type %>% 
  purrr::map(
    ~ summarise_mesh10km_typearea_mesh(df_pref01_ksj_l03, .x))
df_plot_h3 <- 
  landuse_type %>% 
  purrr::map(
    ~ summarise_mesh10km_typearea_h3(df_pref01_ksj_l03, .x))
# 3/3 Visualize -----------------------------------------------------------------------
# plot_type(df_plot_mesh[[1]], "sinrin")
# plot_type(df_plot_h3[[1]], "sinrin")
plot_mesh <- 
  purrr::map2(
    .x = df_plot_mesh,
    .y = landuse_type,
    ~ plot_type(.x, .y))
plot_mesh %>% 
  purrr::reduce(`+`) + 
  plot_layout(nrow = 3) +
  plot_annotation(title = "土地利用 面積",
                  subtitle = "3次(1km)メッシュ内の面積を2次(10km)メッシュ単位で集計、可視化",
                  caption = "データ: 国土交通省 国土数値情報 (土地利用3次メッシュデータ 第2.6版 L03-a 平成26年度 http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L03-a.html)\nTwitter: @u_ribo",
                  theme = theme(text = element_text(family = "IPAexGothic")))
ggsave("out_mesh10km.png", last_plot(), width = 10, height = 6)
plot_h3 <- 
  purrr::map2(
    .x = df_plot_h3,
    .y = landuse_type,
    ~ plot_type(.x, .y))
plot_h3 %>% 
  purrr::reduce(`+`) + 
  plot_layout(nrow = 3) +
  plot_annotation(title = "土地利用 面積",
                  subtitle = "3次(1km)メッシュ内の面積を2次(10km)メッシュ化、h3座標での集計と可視化",
                  caption = "データ: 国土交通省 国土数値情報 (土地利用3次メッシュデータ 第2.6版 L03-a 平成26年度 http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L03-a.html)\nTwitter: @u_ribo",
                  theme = theme(text = element_text(family = "IPAexGothic")))
ggsave("out_h3res5.png", last_plot(), width = 10, height = 6)